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ABSTRACT 

We introduce a radiative transfer code module for the magnetohydrodynamical adap¬ 
tive mesh refinement code FLASH 4. It is coupled to an efficient chemical network 
which explicitly tracks the three hydrogen species H,H2,H+ as well as C+ and CO. 
The module is geared towards modeling all relevant thermal feedback processes of 
massive stars, and is able to follow the non-equilibrium time-dependent thermal and 
chemical state of the present-day interstellar medium as well as that of dense molec¬ 
ular clouds. We describe in detail the implementation of all relevant thermal stellar 
feedback mechanisms, i.e. photoelectric, photoionization and H2 dissociation heating 
as well as pumping of molecular hydrogen by UV photons. All included radiative 
feedback processes are extensively tested. We also compare our module to dedicated 
photon-dominated region (PDR) codes and find good agreement in our modeled hy¬ 
drogen species once our radiative transfer solution reaches equilibrium. In addition, 
we show that the implemented radiative feedback physics is insensitive to the spatial 
resolution of the code and show under which conditions it is possible to obtain well- 
converged evolution in time. Finally, we briefly explore the robustness of our scheme 
for treating combined ionizing and non-ionizing radiation. 
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1 INTRODUCTION 

Massive stars continuously shape their environments by 
feedback in the form of radiation, winds and supernovae. 
Not only do they regulate the thermal and velocity struc¬ 
ture of the interstellar medium (ISM) (see e.g. Lopez et al.| 
201]J^ Krumholz^fc Thompson|j2^^ et al.||2013[ 


are technically challenging and important physical processes 
are often omitted or treated in a simplified fashion in order 
to make the calculations computationally tractable. Unfor¬ 
tunately, it appears that if too much of the real physics is 
omitted, then this can significantly influence the effective- 


Walch et al. 2014[ or [Klessen fc Glover||2014| for ageneral 
overview), they also determine the total fraction of their par¬ 
ent molecular cloud that gets converted into stars ([Murray 


ness of the feedback in the simulations (see e.g. Iliev et al. 
2009 Kuiper et al.|[2012|. It is therefore important to con 


20111JDale et al.|2012||Faucher-Giguere et al.|2013l|Hopkins 

et al. 120141). 

However, the impact of the deposited radiative energy 
and momentum on the star formation efficiency, cloud life¬ 
time and turbulent velocity structure is still subject to de¬ 


tinue to search for ways in which to model the physics of 
stellar feedback accurately and efficiently. 

For radiative feedback from massive stars, many differ¬ 
ent approaches have been proposed in the literature, ranging 
from simple approximations (e.g. Dale et al.|2007 Krumholz 
jet al.||2009| jColfn et al^|2013|), to moment-based schemes 


(e.g. Hayes fc Norman|2003 Rosdahl et al.|2013 Skinner & 


bate (compare e.g. the results of Gritschneder et al. 2009 and 

Ostriker||2013 1, ray-tracing techniques (e.g Rijkhorst et al. 

Walch et al.|2012 with those of Dale et al. 12012 1. A powerful 

2006| [Peters et al.[[2010| [Wise & Abel[[2011[) and finally 


approach to explore the interaction of stellar feedback with 
its surroundings are numerical simulations in which we con¬ 
trol the implemented physics and can choose the boundary 
conditions and initial conditions. Simulations of this kind 
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Monte Carlo methods (e.g. Harries 20151 as well as hy¬ 


brid schemes (Paardekooper et al.|2010 Kuiper et al.|2010 1. 


These different approaches have different degrees of sophis¬ 
tication and often wildly different computational costs. In 
general, the higher the accuracy of the method, the more 
costly it is, limiting the applicability of the most accurate 
approaches to situations with only one or a few radiation 
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sources, and simulation durations that cover only a small 
fraction of the total stellar lifetime. 


One important aspect of any model of radiative feed¬ 
back from massive stars is the way in which the radiative 
output of the star is coupled to the thermal and chemical 
evolution of the gas. Although highly sophisticated models 
of the effect of stellar radiation on the chemical composi¬ 
tion and temperature structure of interstellar gas are avail¬ 
able (e.g. Ercolano et al. 2003 Ferland et al. 2013), they 
do not account for the hydrodynamical evolution of the gas, 
and hence cannot be used to study the effect of radiative 
feedback on star formation. At the other extreme, many hy¬ 
drodynamical models of radiative feedback account for the 
chemical state of the gas only to the extent that they dis¬ 
tinguish between ionized and neutral regions. In addition, 
they often adopt a highly simplified treatment of the ther¬ 
mal state of the gas, for instance by assuming that all ionized 
gas is isothermal with a temperature Tion — 10'* K, and all 
neutral gas is isothermal with a temperature Tneut = 10 K 
(see e.g. Dale et al. 2007 Gritschneder et al.||2009 |. Those 
hydrodynamical models that do attempt to move beyond 
this simple approximation and to properly account for the 
coupling between the radiation, the chemical and thermal 
state of the gas and its hydrodynamical evolution are typi¬ 
cally designed to model radiative feedback in high-redshift. 


metal-free gas (see e.g. many of the codes compared in Iliev 


et al. 2006 2009). This gas generally has a low molecular 


fraction, no dust, and relatively simple microphysics in com¬ 
parison to that making up the present-day ISM. 


In this paper we present Fervent, a radiative trans¬ 
fer module designed to model the coupled chemical, ther¬ 
mal and dynamical effects of radiative feedback on the 
present-day ISM. Fervent consists of an adaptive ray¬ 
tracing scheme, inspired by the Moray algorithm of | Wise ^z\ 


Abel (2011), coupled to a simplified chemical network and 


a detailed cooling function. It allows us to self-consistently 
calculate the chemical state of interstellar gas affected by ra¬ 
diation from massive stars, as well as to address the question 
of how the thermal state of the gas is influenced. This is vi¬ 
tal for an accurate modelling of Hll region dynamics, as it is 
the thermal structure of the gas that determines the strength 
of the hydrodynamical response to the deposited radiative 
energy, which in the end constrains the effectiveness of the 
feedback. In other words, our scheme has the advantage that 
it connects the underlying chemical processes directly to the 
hydrodynamics without an intermediate model for the ther¬ 
mal structure of the gas. Ultimately, we aim to quantify the 
range and strength of stellar radiative feedback from the 
ground up, by including all of the important chemical pro¬ 
cesses that influence the behaviour of the gas. 


Our paper is structured in the following way. In Sec¬ 
tion we present the physical details of the model and an 
outline of our ray-tracing scheme. Tests of our method are 
described in section and we summarize our findings in sec- 
tion|^ Additional technical details regarding the ray-tracing 
scheme are given in Appendix [A] 


2 RADIATIVE TRANSFER PHYSICS AND 
NUMERICS 


The rate of radiative energy transfer dE/dt in erg s ^ at 
every position in space x and time t is given by 


dE = cos On dlln dAn du dt (1) 


where I{y, n, x, t) is the specific intensity in 
erg s“^ Hz“^ cm“^ sr“^ of the radiation field seen in 
the direction of n in the frequency range v to v + du. In 
addition, the amount of locally transferred energy depends 
on how much of the emitting source area An is seen in the 
solid angle fin. 

The evolution of I is modeled by a transport equation 
of the form: 

-I- n • VI = -K(i/,x,t) /-I- j(i/,x,t) (2) 

The right hand side consists of local source and sink terms, 

i. e. the attenuation of the field by some medium with ab¬ 
sorption coefficient k and spontaneous emission coefficient 

j. On the left hand side, the propagation of the field in time 
and its geometric dilution along the direction n is described. 
It is convenient to reduce equation © to one dimension and 
express the specihc intensity in terms of the photon flux P 
from a single radiation source. By integrating I over the solid 
angle and a closed surface containing the emitting source, 
the photon flux Pv in the frequency range v to v E dv \s 
obtained. 


ld]\ d]\ 

c dt dr 


= -K 


(3) 


where we have neglected the local emission j due to e.g. scat¬ 
tering processes. In practice one cannot always stay in the 
ray-frame, as we will show later for processes other than 
hydrogen-ionizing radiation. 

To allow us to solve for the evolution of the photon 
flux as a function of position and time within our computa¬ 
tional volume, we implement a ray-tracing approach similar 
to that proposed by Wise & Abel (2011) in the magneto- 


hydrodynamical simulation code FLASH 4 (Fryxell et al. 2000 
Dubey et al.|[2008). We briefly describe the method in the 


next section. For implementation details concerning the gen¬ 
eral code structure, parallelization, box-ray intersection cal¬ 
culation and ray propagation, see Appendix]^ Coupling to 
a chemical network explicitly tracking the species CO, C*", 
, H and H 2 is described in section |2.2| and an overview 
of missing physics in the current implementation is given in 
section \Tn\ 


2.1 Ray-tracing 

We follow a split approach to calculate the effects of the ra¬ 
diation held on the gas in our simulation domain. 

First, by tracing rays from all radiation sources outward, we 
gather information on the chemical composition and radia¬ 
tion intensity and use this to calculate the radiative heating 
rate, the atomic hydrogen photoionization rate and the H 2 
photodissociation rate. Fach cell of the mesh has to know 
the gas column and photon hux to and from each source 
illuminating it in order to calculate all quantities necessary 
for input into the chemistry module. In a subsequent step, 
the gas composition and temperature is updated locally on 
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a cell-by-cell basis, with no further input from its surround¬ 
ings. 

Ideally, after a single chemistry update another ray¬ 
tracing step should be performed to check whether the new 
state is converged. However, multiple ray-tracing and chem¬ 
istry updates during a single simulation time step are cur¬ 
rently computationally too expensive. Instead we impose a 
limiter on the size of the hydrodynamical timestep At by 
requiring that the change in the fractional abundance of the 
most radiatively affected species, atomic hydrogen, 

|^3^H,new ^^H,old| ~ (4) 


is less than some factor fu- This implies that the radiation 
field is assumed static during the chemistry step. We use 
the H abundance in our limiter, because the main photo¬ 
chemical processes that we are interested in modelling - the 
photoionization of atomic hydrogen and the photodissoci¬ 
ation of H 2 - both change the H abundance significantly. 
Ensuring that the global evolution time step satisfies 


AU 


^ Atold, 

Axh 


( 5 ) 


where the smallest Atnew of all cells is taken, should there¬ 
fore accurately capture the chemical evolution of the gas for 
small enough fu- Other processes that produce atomic hy¬ 
drogen (e.g. the collisional dissociation of H 2 ) generally do 
not lead to changes in the H abundance that occur rapidly 
enough to trigger the timestep limiter, provided that we 
choose a reasonable value for fu- We discuss appropriate 
values for fu in Section |3.1| After the new thermal and 
chemical state of all cells has been calculated, the hydro- 
dynamical equations are advanced. 

During each simulation timestep, rays originating from 
all radiation sources are generated at the same time. The 
rays are spaced evenly on a spherical surface centered on the 
source locations by utilizing the HEALPix software package 


(Gorski et al. 20051. HEALPix is an equal area decomposi¬ 


tion of a unit sphere with a nested grid structure. The base 
level of refinement consists of twelve pixels, where each pixel 
can be split independently into four nested child pixels. All 
pixels are numbered consecutively which leads to a limit of 
thirteen levels of refinement, where the pixel number just 
hts inside a long integer. This corresponds to a maximum 
angular resolution of 3.2 x 10“^ arcseconds. The indepen¬ 
dent splitting makes it inherently applicable to an adaptive 
mesh refinement (AMR) code such as FLASH 4. 

After generating the initial sphere of rays, intersections 
between them and the cubic cells of the simulation mesh 
are calculated during the ray-tracing step. Rays are repre¬ 
sented by an equation of a line in vector form p(r) = s + rn, 
with the unit direction vector n, source location s and total 
traversed distance r. Each time a ray enters a new cell the 
intersection distance dr from the entry to the exit point is 
calculated. To compute the radiative heating rate and the 
various photochemical rates, we need to know for each ray 
and each cell the value of the ray segment dr in that cell, as 
well as the total column densities of atomic and molecular 
hydrogen, Ah and Nu 2 , traversed by the ray between the 
source and the cell. These column densities are given by 


k 

A^h,H 2 = riidri, (6) 

i=0 


where i = 0 represents the cell containing the point source, 
k is the current cell, dri is the ray segment in cell i, a:H,H 2 (*) 
are the fractional abundances of H and H 2 , respectively, in 
cell i, and rii is the number density of H nuclei in cell i- These 
quantities are computed as we move along the rays and are 
saved in the mesh data structure in the same way as any 
other hydrodynamical state variable. For the propagation 
of the radiation field, the gas columns are carried to the 
next neighboring cell as a ray property (see Appendix [A| for 
details). 

For a full sampling of the radiation field, each mesh cell 
has to be intersected by at least one ray. At this point, we 
make use of the HEALPix software’s ability to split pixels by 
assigning each ray a HEALpix level I and index h- We use the 
splitting criterion introduced in Wise & Abel (20111, which 
calculates the ratio of the mesh cell size and the pixel size 
at a distance r, 


f _ ^ _ d? 

where A is the face area of a cubic cell with side lengths 
d, D is the solid angle associated with a single pixel at the 
current refinement level and Apixei = 12 x 2^ x 2* is the total 
number of pixels on refinement level L Provided that fn > 1, 
this criterion guarantees that each cell is sampled by at least 
one ray. Conveniently, this criterion also accounts for regions 
refined during the course of the simulation, where instead of 
a global uniform grid with constant area A, local regions are 
resolved by a finer mesh. This results in smaller cell sizes and 
rays are split to adapt to this difference. Multiple rays per 
cell are treated by summing up all individual contributions 
to the photochemical rates and heating rates, i.e. 


^ rays 
^tot — ^ ^ 

i=0 


( 8 ) 


without any additional weighting. In the same vein, multiple 
sources are taken into account. We do not distinguish be¬ 
tween rays from different sources, and so the order in which 
the individual contributions are calculated is not controlled. 

This ray-tracing approach is intrinsically serial as new 
rays have to be generated during the traversal of the sim¬ 
ulation domain. The total number of rays is therefore un¬ 
predictable at the beginning of a simulation step. In con¬ 
trast, the most efficient parallelizing schemes expect that all 
ray information is available from the beginning of the ray¬ 
tracing step and parallelize ray-cell intersection calculations, 
see e.g. Rijkhorst et al. (20061. However, this has the disad¬ 


vantage that for a large number of sources, all rays from all 
sources have to be generated at the same time, consuming 
large amounts of memory. In addition, this introduces an 
overhead in communication as rays have to be known glob¬ 
ally everywhere to account for possible local contributions. 
A simple way to reduce the required amount of communica¬ 
tion is a pruning step where only the rays that actually enter 
the local domain are to be communicated. Another way to 
prepare the rays could be to pre-generate all intersections in 
a tree-walk step. 

In our scheme each domain containing a source starts 
to locally ray-trace. If all sources are distributed evenly over 
all domains this approach is efficient, but in most astrophys- 
ical applications sources are expected to cluster. However, 
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as rays are traced from the source outward, once the op¬ 
tical depth T becomes too large they can be terminated, 
essentially introducing an on-the-fly pruning. This greatly 
reduces the cost of the ray-tracing step. 


2.2 Photochemistry and radiative heating 


In this section, we describe how we account for the effects 
of the radiation field on the chemistry and thermal balance 
of the gas. The processes that we currently treat are the 
photoelectric heating of the gas by the dust, the photodisso¬ 
ciation of the molecules H 2 and CO, the pumping of excited 
rotational and vibrational states of H 2 by ultraviolet (UV) 
photons, and the photoionization of atomic and molecular 
hydrogen. To do this, we couple the reaction and heating 
rates gathered during ray-tracing to the NL97 chemical net¬ 


work described in Glover & Clark 1 2012a I; note that details 


of the implementation of this network in FLASH can be found 
inf 


Walch et al. (20141. 


Our main concern is radiation from stars, which for sim¬ 
plicity we approximate by a blackbody spectrum i.e. 
we neglect any absorption or emission line features. We di¬ 
vide this spectrum into four energy bins: 5.6 eV < Ess < 

11.2 eV, 11.2 eV < £;ii .2 < 13.6 eV, 13.6 eV < Bias < 

15.2 eV and E 15 . 2 + > 15.2 eV. Using a larger number of 
energy bins would allow us to account in more detail for 
the frequency-dependent opacity of the gas, but at the cost 
of slowing down the entire calculation. We have found that 
using four bins appears to give us the best trade-off be¬ 
tween speed and the accuracy with which we can model the 
thermal balance of the gas and the details of the hydrogen 
chemistry. For the tracked species we explicitly calculate the 
rates of the reactions: 


( 9 ) 


where the subscript number denotes the lower end of the 
photon energy bin and the plus sign that the bin extends to 
infinity. Reaction (i) describes the dissociation of H 2 through 
excitation by far ultraviolet photons. We assume that the 
only photons responsible for this process are those in the 
energy range 11.2 < hu < 13.6 eV. Lower energy photons 
cannot dissociate H 2 molecules that are in the vibrational 
ground state and hence become significant only if we are con¬ 
cerned with the photodissociation of H 2 in hot, high density 


i) 

H 2 

-I- 711.2 

^ H 

-f H 

ii) 

H 

+ 713.6 

^ H+ 

- 1 - e“ 

iii) 

H 

+ 715 . 2 + 

-+ H+ 

- 1 - e“ 

iv) 

H 2 

+ 715 . 2 + 

^H+ 

- 1 - e“ 


Glover 


20151, condi- 


gas (T > 2000 K, n > lO'^ cm“'^; see 
tions in which collisional dissociation of H 2 is in any event 
likely to dominate. Photons with energies hiy > 13.6 eV can 
dissociate H 2 , but have a much larger probability of being 
absorbed by atomic hydrogen resulting in photoionization. 

Reactions (ii) and (iii) account for the photoionization 
of atomic hydrogen by the radiation field. For reasons which 
will become clear below, we distinguish between photons 
with energies 13.6 < hu < 15.2 eV that can ionize atomic 
hydrogen but not molecular hydrogen, and photons with en¬ 
ergies hu > 15.2 eV that can ionize both H and H 2 . Finally, 
reaction (iv) accounts for the photoionization of H 2 . This 
produces ions, but this chemical species is not directly 
tracked in the NL97 chemical network. Instead, we assume 


that all of the Hj ions are destroyed by dissociative recom¬ 
bination 

-f e"H-h H, (10) 

so that a single photoionization of H 2 ultimately results in 
the production of two hydrogen atoms. Since the rate coef¬ 
ficient for this reaction is large, and we expect there to be 
a high electron density close to the ionization front, it is a 
good approximation to treat this process as effectively in¬ 
stantaneous (i.e. we do not expect to build up a significant 
abundance of Hj). 

We note that if one does not account for the effects 
of photoionizing radiation on the H 2 then it is easy to ob¬ 
tain unphysical results in highly molecular gas, since in this 
case the bulk of photons with energies larger than 13.6 eV 
can freely escape from a fully molecular medium, where the 
abundance of atomic hydrogen and thus the effective optical 
depth are zero. In that sense, our treatment here also forms 
a numerical closure. We considered several alternative ways 
of dealing with this problem, such as slowly ramping up the 
strength of the radiation sources (allowing a photodissoci¬ 
ation region (PDR) time to form before the gas starts to 
become ionized), or artificially creating an initial PDR, but 
we found that the behaviour of these alternatives depended 
too strongly on the spatial resolution of the simulation and 
was hard to parameterize. 


2.2.1 Ionizing radiation: chemical effects 

Photons with an energy E,, = hu greater than 13.6 eV can 
ionize atomic hydrogen. A star with a blackbody spectrum 
and effective temperature of Teff emits 

Pion = / 47 rR* --- du= Au ( 11 ) 

photons capable of ionization per second, where h is the 
Planck constant and R* is the radius of the star. We assume 
the photon flux to be constant over one timestep, which is 
reasonable since on the scales that we are interested in mod¬ 
elling, the hydrodynamical timestep is orders of magnitude 
smaller than the lifetime of even the most massive O-type 
star. In this case, equation ||^ reduces to 

^ = -K Pan (12) 

dr 

where we insert the ionizing photon flux Pion for P,^. At 
this point two approaches are possible to calculate the at¬ 
tenuated flux in a given cell. One option is to use the total 
distance r from the source location at r = 0 to the current 
position 

r ^ = r -K dr', (13) 

J Pq ion J 0 

from which it follows that 


P(r) = Po e-", 


(14) 


where Po = Pion(r = 0), P = Pion(r), and r = k dr' is 
the total optical depth along the ray. Alternatively, we can 
use the ray segments dr. 


L 


P+dP 


- dP- = 

p, ion 


/ 


r+dr 

—K dr' 


(15) 













Fervent: Chemistry-coupled radiative transfer 5 


to arrive at 


dP = P 


(16) 


which gives the change in P over a single cell. Each treat¬ 
ment of the photon flux requires different properties to be 
carried along the ray, either Pq or the incrementally atten¬ 
uated flux Pi+i — Pi — dPi, where i denotes the current cell 
index. Depending on how photons of the radiation field are 
absorbed by the medium, the first or second formalism is 
more advantageous. If there is a reaction that maps a single 
photon to a transition from one species to another and all 
attenuated photons lead to this transition, dP provides an 
easy prescription for equating all absorbed photons to the 
change in the appropriate species. In this case the scheme is 
photon conservative and the incremental change in optical 
depth dr = crrisdr can be expressed in terms of a reaction 
cross-section a, the number density of the affected species 
Us and the ray segment dr. 

In comparison, if a transition to a new species only oc¬ 
curs for a fraction of the absorbed photons, which is deter¬ 
mined by some non-local quantity such as the total column 
of gas to the current cell position, dP cannot be photon 
conservative. In that case the transition rate has to be de¬ 
termined from the flux at the source Pq. 

To model ionizing radiation, we follow [Wise fc Abel| 


120111, who use the second approach and attenuate the pho¬ 


ton flux based on ray segments dr. We start by converting 
Po to the absolute number of ionizing photons emitted in a 
timestep At at the source position 


Nion{r = 0) = Po At. 


(17) 


to determine how many photons take part in each process. 
One approach would be to simply divide the photons into 
two groups by using the relative abundances of Ho and H. 
However, this approach fails to account for the fact that the 
photoionization cross-section of Ho, (thj, is different from 
that of H, cth. We therefore divide up the photons in a way 
that accounts for both the relative abundances of Ho and 
H and the relative size of their ionization cross-sections. We 
denote the number of photons in the P 15 . 0 - 1 - energy bin that 
are removed from the ray in a given cell by Ho photoion¬ 
ization as and the number that are removed by H 

photoionization as These are given by 


diV'®" 


dAi5.o+ 

dAi5.o+ 


/r^Ho 

(^Ho/r -I- Xh) ’ 
Xh 

{xH2fr -I- a:H)' 


( 20 ) 


Here, dAi 5 . 2 -i- is the total number of photons in the P 15 . 2 + 
energy bin that are absorbed in the cell in question and 

f — (^^ 2 ) , 21 ') 

where {cyi^) and are intensity-weighted mean cross- 

sections for H and Ho ionization by photons with E > 
15.2 eV. These are given by the expressions 


(CTH2(Piff)) = 


15 . 2 + 


(^eff)) = 


r 

J 

f 

J 1/1 


2 

anOBu 

hv 


df 

J ui 

' J 


B, 

hu 


dv. 


B. 

hu 


( 22 ) 


du, 


This number of photons is then divided evenly amongst 
A^pixei initial rays, so that each ray starts with Pq At/Apix 
ionizing photons associated with it. This initial set of rays 
is set up so that each ray has a direction vector centered 
on (and normal to) the associated HEALPix pixel, pointing 
away from the source. We then walk along each ray. As we 
enter each new cell along a given ray, we reduce the num¬ 
ber of ionizing photons currently associated with that ray, 
Nion{r), by an amount dWon that is given by 

dWon = Aio„(r) (1 - e-'^"). (18) 


We next exploit the fact that dWon is a dimensionless quan¬ 
tity to allow us to formulate the ionization rate in a given 
cell in a way that is independent of spatial resolution. This 
is done by calculating the total number of neutral hydrogen 
atoms in a cell, using the local number density of H atoms, 
riH, and the volume of the cell, Keii = cT^, where d is the 
length of one side of the cell and m is the dimensionality of 
the mesh. The ionization rate then follows as 


^ion — 


dWc 


w-nWeii At 


(19) 


This formulation of the ionization rate ensures, by construc¬ 
tion, that the number of ionizing photons absorbed is equal 
to the ionization rate per unit volume integrated over the 
volume of the cell, i.e. it is photon conservative. 

For photons in the P 13.6 energy bin, this is all we need 
to worry about, at least as far as the photoionization rate is 
concerned. For photons in the Pi 5 , 2 + energy bin, however, 
things are more complicated, since these photons can ion¬ 
ize either atomic or molecular hydrogen. We therefore have 


where Peff is the effective temperature of the stellar source, 
Bi/ = BviTgff) is the corresponding Planck function, auO 
and (JH 2 O are the frequency-dependent ionization cross- 
sections for H and H 2 and hui 5.2 = 15.2 eV. For a^iu), we 
use the simple approximation ( Osterbrock|1989 1 


I , / 1 ^ 13.6 

aniy) = (Jufi j , 


(23) 


as this is sufficiently accurate for the photon energies of in¬ 
terest here. For aii 2 {Cj a piecewise constant cross- 

section to the analytical results of |Liu fc Shemansky (20121 
for energies in the range 15.2 < E < 18.1 eV, obtain¬ 
ing the values shown in Table At higher energies, we 
assume that the H 2 ionization cross-section falls off as 
= o-H 2 ,high(ffi 8 .io/j^)^, where hi^ig.io = 18.10 eV and 
o'H2,high = 9.75 X 10“^® cm^. At these energies, our adopted 
values differ from those suggested by experiment by no more 
than 20% ( [Chung et al.|[T9^ |. The behavior of au{u) and 
crH 2 (i^) as a function of energy is illustrated in the left-hand 
panel of Figure 

Once we have computed dA^?'^ and dAjQ^'^, it is 
straightforward to then calculate the rate at which H 2 is 
destroyed by photons in the Bi 5 , 2 + energy bin: 


1 . 15.2 _ 


n-HaKell At 


(24) 


Similarly, the rate at which atomic hydrogen is photoionized 
by photons in the B 15 . 2 + energy bin is simply 


, 15.2 _ dNtfC 

nnHceiiAf 


( 25 ) 
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Figure 1. Left: frequency dependent cross-sections a for the ionization of H and H 2 by photons with energy Ep^,. The hydrogen 
photoionization cross-section is denoted by a dotted line up to an energy of 15.2 eV, the H 2 photoionization threshold, and by a 
dashed line at higher energies. The solid line denotes the H 2 photoionization cross-section Right: intensity-weighted mean H and H 2 
photoionization cross-sections for a range of different effective temperatures. The dotted line shows the intensity-weighted mean 

photoionization cross-section for H in the i?i 3.6 energy bin, and the dashed line shows the same quantity in the Ei 5 , 2 -|- energy bin. 

The solid line shows the intensity-weighted mean photoionization cross-section for H 2 in the Ei^, 2 + energy bin; the corresponding 

value in the lower energy bin is zero. 


Table 1. Piecewise Cross-Section for Ionization of H 2 


a (Mb) 

Emin (eV) 

^max (sV) 

0.00 

13.60 

15.20 

0.09 

15.20 

15.45 

1.15 

15.70 

15.95 

3.00 

15.95 

16.20 

5.00 

16.20 

16.40 

6.75 

16.40 

16.65 

8.00 

16.65 

16.85 

9.00 

16.85 

17.00 

9.50 

17.00 

17.20 

9.80 

17.20 

17.65 

10.10 

17.65 

18.10 


18.10 

00 


Note that with this formulation of the H 2 and H photoioniza¬ 
tion rates, we guarantee that the total number of photoion¬ 
izations balances the total number of photons absorbed in 
the cell, regardless of the relative sizes of nu and nua. 

In order to compute the total photoionization rate of 
atomic hydrogen, we need to account for the ionizing pho¬ 
tons in the ifis.e energy bin (reaction ii) as well as those 
in the Ei 5 , 2 + energy bin (reaction iii). To compute the rate 
due to the softer photons, which do not couple to H 2 , we 
simply compute 


. 13.6 _ dAfl3.6 

“nHl/cellAf 


(26) 


where dAfi 3,6 is the total number of photons in the Eis.g 
energy bin that are absorbed in the cell. The total pho¬ 
toionization rate then follows trivially: 

k- — -I- (27) 

'Tion — '^lon I ^lon * ‘ / 


Another aspect of the problem that needs to be treated 
with care is the fact that if there is no atomic hydrogen 
in the cell, then none of the photons in the ifis.e energy 
bin can be absorbed. If we therefore compute the number 
of photoionizations from this bin hrst and from the Tlis. 2 - 1 - 


bin second, we arrive at an unphysical scenario: the lower 
energy photons merely propagate through the cell without 
any being absorbed (because at that point in the calculation 
there is no H present) and the higher energy photons then 
destroy some of the H 2 , creating some atomic hydrogen, but 
too late for this to affect the lower energy ones. To avoid this 
problem, we simply ensure that we account for the effects 
of H 2 dissociation before dealing with the Eis.e energy bin. 
If the H 2 formation time is long compared to the H 2 disso¬ 
ciation time (which is almost always a good approximation 
in the clouds that we model), then the change in the atomic 
hydrogen number density due to the destruction of H 2 by 
energetic photons can be written as 

diV).®'"+ 

dnu = 2fcdisnH2 At = 2——^5—. (28) 

hcell 

When computing we therefore do not use nn, the num¬ 
ber density of atomic hydrogen at the start of the timestep, 
but instead use the improved estimate n'u = uh + dnn that 
accounts for the destruction of H 2 by energetic photons dur¬ 
ing the timestep. If nn is very small or zero, then this pro¬ 
cedure allows us to avoid the coupling problem described 
above. On the other hand, if nn is larger (i.e. if most of the 
H 2 in the cell has already been destroyed), then it becomes 
only a minor correction. 

To complete our specification of how we determine fcion, 
we need to describe how we calculate dAii 5 , 2 -i- and dN;),^'®. 
We start by writing them as 

dNi3,6 = Nis.e [1 - exp(-dri3,6)], (29) 

and 

dNl5.2+ = ^15.2-1- [1 — exp( —dTi5,2+)] , (30) 

where Nis.e and Aii 5 . 2 -i- are the numbers of photons in the 
Ei 3 ,e and E 1 S. 2 + energy bins at the point where the ray 
enters the grid cell, and dris.e and dTi 5 . 2 + are the change in 
the optical depth of the two energy bins taking place across 
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the cell. These are given by 


and 


dri3.6 = (CTH^'®)nHdr, 


dri5.2+ = [(crH®'^''")nH + (crH 2 )nH 2 ] dr, 


where 


(alfHnC) = r 

J 


(TuOBl^ 

hv 


dv / 


A 


(31) 


(32) 


(33) 


is the intensity-weighted mean photoionization cross-section 
for H in the i?i 3.6 energy bin, and where and {(TH 2 ) 

have the same valnes as before. When computing {uh®'®), 
we again use the simple approximation for the frequency- 
dependent atomic hydrogen photonionization cross-section 
given in equation ( [^ . 

Finally, once we have computed fcion for each ray passing 
through the grid cell, we sum the contributions from all of 
the rays to get the total photoionization rate: 




^ ) ^ion. 


(34) 


This is then passed to the chemistry module, which solves 
the following rate equation for the H’*' abundance: 


dXH+ 

dt 


— CcineXn (C*xr -f tJcr -f ^ion,tot) :rH C*cr ^H2 

CTB Ue3:]3+ • 


(35) 

Here, in addition to photoionization, we also account for 
the ionization of atomic hydrogen by cosmic rays (Ccr), X- 
rays (Cxr), and its collisional ionization by electrons (Cd). 
Additional ions are also produced by the dissociative 
ionization of H 2 by cosmic rays (Ccr), while ions are 
removed from the gas by case B radiative recombination 
(qb) and recombination on the surfaces of dust grains (hgr). 
Full details of our implementation of all of these processes 
can be found in Glover et al. ( |2010 1 and Walch et al. 120141, 
and so for brevity we do not repeat them here. Finally, we 
note that when computing the electron abundance, we do 
not simply account for the electrons coming from ionized 
hydrogen, but also those from atomic carbon, i.e. Xe = 2:h+ + 
Xc+- 


2.2.2 Ionizing radiation: thermal effects 


Each time a photon unbinds an H 2 molecule or ionizes a 
hydrogen atom, its excess energy is deposited in the sur¬ 
rounding gas in the form of heat. The average deposited 
energy (E) is given by the expression 


{E) 






a,, dv 


(36) 


where the appropriate frequency-dependent cross-sections 
<^H 2 ,i^ and lower bounds should be inserted to 
yield the mean energies deposited by photoionization of H 
due to photons in the i?i 5 . 2 -i- energy bin, photoion¬ 

ization of H due to photons in the E 13.6 energy bin, {E^Ciy 
and the dissociation of H 2 due to photoionization, {Edia}. 



Tft (K) 


Figure 2. Top: average excess energy per H ionization, for ion¬ 
izing photons in the energy bins Fis.e (dotted) and F 15 . 2 + 
(dashed), plotted as a function of the effective temperature of 
the stellar source. The solid line shows the average excess energy 
per H 2 ionization. Bottom: Ratio of the heating rate per H ion¬ 
ization to that for H 2 ionization, considering only photons in the 
F 15 . 2 + energy bin. 


Figure shows how these quantities vary as a function of 

the effective temperature of the source. 

Given the average energy per ionization, it is then 

straightforward to compute the corresponding heating rate 

per unit volume. We have 

Ti _ r T 13.6 / pis.6\ I I 15.2 / pl5.2\ 1 

-L ion — [^ion \-^ion / “r ^ion \-^ion /J (37) 

Tdis ~ [^dis (-E'dis)]'^H 2 • 

These heating rates are accounted for when we determine the 
change in the internal energy of the gas during the timestep, 
as we describe in Section [2.2.61 below. 


2.2.3 Non-ionizing radiation: chemical effects 


Even though a minimum photon energy of about 15.2 eV 


is needed to ionize molecular hydrogen (Liu & Shemansky 


2012 1, less energetic photons are able to dissociate H 2 by a 
two-step process ( [Stecher &: Williams|1967| . A photon with 
energy larger than 11.2 eV is able to electronically excite 
H 2 to the or C states (also known as the Ly¬ 

man and Werner states). As there are a large number of 
different bound rotational and vibrational levels in both the 
electronic ground state and the Lyman and Werner excited 
states, these electronic transitions occur via a series of dis¬ 
crete lines that together make up the Lyman and Werner 
band systems]^ In most cases, the excited H 2 molecule de¬ 
cays back into a bound rotational and vibrational level in 
the electronic ground state. However, roughly 15% of the 


^ Since these sets of lines overlap in terms of frequency, it is 
common to refer to the range of energies corresponding to the 
combined set of lines as the Lyman-Werner band. 
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time, the decay occurs instead into the vibrational contin¬ 


uum, resulting in the dissociation of the molecule | Draine & 
Bertoldijl996 1 . 


In contrast to photoionization, which is a continuum 
process, H 2 photodissociation is inherently a line-driven pro¬ 
cess. This means that it is not possible to treat it in the same 
way that we treat photoionization, as it is important to ac¬ 
count for the effects of H 2 self-shielding. Consider a beam of 
radiation propagating along a ray that passes through some 
molecular hydrogen. As our beam of radiation propagates 
through the gas, the photons that are most likely to be ab¬ 
sorbed by the H 2 are those with frequencies corresponding 
to the center of one of the strong Lyman-Werner absorp¬ 
tion lines. As these photons are removed from the beam, 
the probability of additional photons being absorbed by H 2 
decreases and so the photodissociation rate declines. This ef¬ 
fect is known as self-shielding, and it is frequently the domi¬ 
nant process responsible for shielding molecular gas from the 
effects of photodissociation. Importantly, the rate at which 
the Lyman-Werner absorption lines become optically thick 
depends on the strength of the lines: we can easily find our¬ 
selves in a situation in which the strongest lines are optically 
thick, while the weaker lines remain optically thin. In addi¬ 
tion, even once the strongest lines become optically thick 
at their centers, the wings of the line often remain optically 
thin. As a result, the photodissociation rate at a given point 
along the ray depends in a non-trivial way on the H 2 col¬ 
umn density and the velocity structure between that point 
and the source of the radiation. Consequently, we cannot 
easily use the ray segment method to treat the effects of 
photodissociation, as the reduction in the photodissociation 
rate that occurs as our ray crosses any given grid cell does 
not depend simply on the column density of H 2 within that 
cell, but also on the total column density of H 2 between the 
cell and the source. 

We therefore employ a modified version of equation ( |13[ ) 
to treat the H 2 photodissociation rate: 


fcuv = iVii. 2 (r = 0)e-"<i^5^^. (38) 


Here Aii. 2 (r = 0) = Pii. 2 {r = 0)At is the absolute number 
of UV photons with energies 11.2 < E < 13.6 eV emitted by 
our source in a timestep At, auv is the effective photodis¬ 
sociation cross-section, and Td is the mean optical depth in 
the Lyman-Werner band due to dust absorption. 

The effective photodissociation cross-section, ctuv, is 
defined as the ratio of the photodissociation rate D to the 
photon flux in the Lyman-Werner bands, F, i.e. cruv = 
DjF. We evaluate this using the photodissociation rate in 


the optically thin limit, D = 5.18 x 10 s ^ 


_ [Rollig et al. 

|2007[ ), which was derived for the |Draine| ( |1978P interstellar 
radiation field (ISRF). The ISRF approximates the aver¬ 
age spectrum found in the interstellar medium (ISM) orig¬ 
inating from reprocessed starlight. With the photon flux in 
the Lyman-Werner bands that o ne has for the same ISRF , 
namely F = 2.1 x 10^ s“^cm“^ ( Draine fc Bertoldi]|l996 1 , 
we find cruv, thin = 2.47 x 10~^® cm'^ in optically thin gas. 
As the H 2 column density increases, however, the gas starts 
to become optically thick in the Lyman-Werner lines. This 
reduces D but has little effect on F until we reach extremely 
large H 2 column densities. To account for this, we write the 


effective photodissociation cross-section as 


cruv — CruV.thin/shd, 


( 39 ) 


where /ahd is a self-shielding function that parameterizes 
how D decreases as the H 2 column density increases. For this 
shielding function, we adopt the expression given in |Drain^ 
|fc Bertoldi| ( |1996l ): 


_ 0.965 0.035 

(1 -I- xjb^Y "*"(! + 

exp 1^—8.5 X 10~'’(1 -|- . 


(40) 


Here, x = Ah2/(5 x 10^"^ cm~^), where Ana is the column 
density of H 2 along the ray between the cell and the source, 
and 65 = 6 / 10 ® cms“^, and 6 parameterizes the change in 
shielding behavior for lines broadened due to the Doppler 
effect. For simplicity, in the current version of Fervent we 
do not relate 6 to the temperature or velocity distribution 
of the H 2 along the ray, but instead simply set 65 = 1 . 0 . 
We plan to relax this assumption in future work, but do not 
expect that it will make a large difference to the outcome of 
our models. 

We assume that the dust has properties characteristic 
of the diffuse ISM. In this case, Td = qAv with 7 = 3.5, and 
the visual extinction Av is related to the total hydrogen 
column density via ( [Draine k. Bertoldiyi996[ ) 


— + 2Ah 2 + Aue r 

1.87 X 1021 cm-2 


(41) 


where /d/g is the dust-to-gas ratio, normalized to the 
value in the local ISM. The effective dust absorption cross- 
section per hydrogen nucleus is therefore Ud = rd/(1.87 x 
10 ^^) cm 2 _ ^ ^ cm 2 , jg much smaller 

than (Tuv.thin, and dust attenuation only becomes impor¬ 
tant when the H 2 fraction is low or when there is a very 
large column density of H 2 between the cell and the source. 
We are therefore justified in treating the effects of H 2 self¬ 
shielding prior to those of dust absorption when accounting 
for the photons absorbed in the grid cell. Although not im¬ 
portant in the present context, this becomes important once 
we consider photoelectric heating, as we have to take care to 
avoid the double-counting of photons: those absorbed by H 2 
are not available to be absorbed by dust and hence cannot 
contribute to the photoelectric heating rate. We return to 
this point in Section [2. 2. 5| below. 

Putting everything together, our final rate equation for 
the molecular hydrogen abundance reads: 


da:H 2 

dt 


— (C'cl.H n-H + Ccl,e Ue -|- Cd.Ha «H2 ) a;H2 

— (fccr + fciSRF + kuv + kdis) X 1 I 2 + RdXn, 


(42) 


where Cci,h, Cci,e and Cci.h2 are the rate coefficients for 
the collisional dissociation of H 2 by H atoms, electrons and 
H 2 molecules, respectively, ka is the rate at which H 2 is 
destroyed by cosmic rays, fcisRF is the rate at which it is 
photodissociated by the ISRF, and Rd is the rate at which 
it forms on dust grains. For more details on all of these 
processes and their implementation in FLASH, we refer the 
reader to [Glover et al.j ( [2010[ ) and [Walch et al.[ ( [2014[ ) . 

Finally, we should also briefly consider the other main 
chemical effect of the non-ionizing photons in our model, 
namely the photodissociation of CO. The CO photodisso¬ 
ciation threshold is 11.09 eV, similar to that of H 2 , and so 



































Fervent: Chemistry-coupled radiative transfer 9 


in principle both CO and H 2 compete for the same set of 
photons in the -Bn .2 energy bin. However, the details of the 
CO self-shielding are quite different from those of the H 2 
self-shielding, and in principle one should also account for 


Following the notation of Draine & Bertoldi (19961, the con- 


the shielding of CO by H 2 and vice versa (see e.g. Visser 
et al.|2009l. Properly accounting for these effects would add 


significant complexity to the code, but is of questionable 
utility given our highly simplified treatment of CO chem¬ 
istry, which is known to overproduce CO by a factor of a 
few ( Glover Sz Clark]|2012a l. Therefore, for the time being 
we have made the major simplification of assuming that the 
CO photodissociation rate scales directly with the total H 2 
dissociation rate, i.e. 

fcco = /conv(A:uv + kdis). (43) 

The conversion factor is taken to be /conv = 3.86 (Rollig 


|et al.|[2007[ ), the ratio of the CO to H 2 photodissociation 
rates in optically thin gas. We note also that it is important 
to include the fedis term to ensure that CO cannot survive in 
highly ionized regions, as this behavior would be unphysical. 
This approximation produces reasonable results in highly ir¬ 
radiated gas, but predicts that fcco falls off too rapidly in the 
Ay ~ 0.1-1 regime. As our tests in Section [3.4| demonstrate, 
this can result in an order of magnitude overestimation of 
the total CO column density in optically thick gas. Our cur¬ 
rent scheme is therefore suitable for predicting whether a 
given molecular cloud or clump is likely to be CO-bright or 
CO-faint, but not for making detailed predictions of e.g. CO 
or C”*" emissivities. Note, however, that our overestimation of 
the CO abundance is unlikely to significantly affect the ther¬ 
mal structure of the gas. CO is a more effective coolant than 
C^ only in very cold gas (see e.g. Glover & Clark 2012b I, but 
the region where we overestimate its abundance is warmed 
by the same radiation responsible for destroying the H 2 and 
CO (see e.g. Section [3^. 


2 . 2.4 Non-ionizing radiation: thermal effects 

Each Lyman-Werner band photodissociation of an H 2 
molecule deposits (Buv) of energy as heat. Typically, we 
hnd that (Buy) — 0.4 eV ( Black &: Dalgarno|1977 |. In addi¬ 
tion to this, UV photons can also heat the gas by indirectly 


exciting the vibrational levels of the H 2 molecule (Burton 
let al.|1990| . As we have already mentioned, the absorption 
of a Lyman-Werner band photon results in photodissocia¬ 
tion only around 15% of the time. The rest of the time, the 
electronically excited H 2 molecule decays back to a bound 
ro-vibrational level in the electronic ground state. A small 
fraction of these decays put the molecule back in the u = 0 
vibrational ground state, but in most cases, the H 2 molecule 
is left with a considerable residual internal energy in the 
form of vibrational excitation. In low density gas, this en¬ 
ergy is simply radiated away as the H 2 molecule undergoes 
a series of radiative transitions that eventually place it back 
in the vibrational ground state. In dense gas, however, col- 
lisional de-excitation can be more effective than radiative 
de-excitation, and in this case most of this energy is redis¬ 
tributed as heat. 

The rate at which vibrationally-excited H 2 is produced 
- often referred to as the UV pumping rate - is related to 
the H 2 photodissociation rate by 


version factor between the two rates is given by 


/p 


1.0 (Pdiss) (Pret) 
(Pdiss) 


(45) 


where (pdiss) is the mean photodissociation probability and 
(piet) is the mean probability that the molecule decays di¬ 
rectly back into the u = 0 level of the electronic ground 
state. These values vary a little depending on the shape of 
the incident spectrum, the degree of H 2 self-shielding that is 
occurring, and the density and temperature of the gas, but 
in typical ISM conditions, /pump = 6.94 (jDraine & Bertoldij 
T9^ . 

The mean energy that is converted to heat per UV 


pumping event has been calculated by Burton et al. (1990 \ 
and can be written as 


(Bpump) = 2 eV X 


Gdex 


Cd ex Grad 


(46) 


where Gdex is a representative value for the collisional de¬ 


excitation rate, given by Burton et al. (19901 as 


Gdex =10-"^(1.4 e 

1.0 


xnn + 


)Vf 


n, 


(47) 


where T is the temperature in Kelvin, n is the gas num¬ 
ber density in cm“®, and where Grad = 2 x 10“^ s“^ is a 
representative value for the radiative de-excitation rate. In 
warm atomic gas, Gdex ~ Grad for number densities around 
n ~ lO'^ cm~^, while in cold molecular gas an even higher 
density is needed. It is therefore clear that this process is 
important only in relatively dense PDRs. 

Finally, putting these two contributions together, we 
can write the total heating rate per unit volume due to the 
absorption of Lyman-Werner band photons as 

Fuv = (fc pump (Bpump) + fcuv(Buv)) nHa • (48) 


2.2.5 Photoelectric heating 

The final important effect that we need to account for is 
photoelectric heating. Photons with energies greater than 
5.6 eV are able to dislodge electrons from dust grains, and 
these electrons go on to collisionally heat the surrounding 
gas. We assume that the only photons that contribute to the 
photoelectric heating rate are those in the Bs.e and Bn .2 en¬ 
ergy bins, since all of the photons in the two higher energy 
bins will be consumed in the ionization of H or (if energetic 
enough) H 2 . The lower end of the photo-electric heating en¬ 
ergy range is not well defined and we adopt for simplicity 


the cut-off used in Draine & Bertoldi (19961 for the ISRF. 


We follow the prescription by Bakes & Tielens (19941, 


who assume a dust grain size distribution that extends down 
to structures as small as polycyclic aromatic hydrocarbons 
(PAHs) and calculate the heating rate Fpe in erg cm“® s“^, 
which takes into account electron and ion recombination 
cooling 


fcp 


— /pumpfcuV- 


(44) 


Fpe = 1.3 X 10“ 


^Gq gTl. 


(49) 










































10 C. Baczynski, S. C. O. Glover, R. S. Klessen 


The photoelectric heating efficiency is given by 


4.9 X 10” 


1 + 4 X 10-3(Goyr/ne<^pah)°'^3 
3.7 X 10"^(r/10^)° '^ 

1 + 2 X 10-4(Goyr/ne</<pah)’ 


(50) 


where Go is the strength of the local, incident interstellar 
radiation field normalized to the integrated Habing held 


(Habing 19681, T is again the temperature in Kelvin and 
rie is the electron number density in units of cm“®. The 


version of t shown here is taken from Wolhre et al. (20031, 


and includes a dimensionless scaling parameter ())pah that 
was not included in the original Bakes & Tielens (19941 


prescription, as well as a 30% larger pre-factor in equation 


(49l. These modihcations were introduced by Wolhre et al. 


to account for the fact that observations of the ISM by the 
Infrared Space Observatory imply the presence of a larger 
abundance of PAHs than was assumed in (Bakes fc Tielensl s 
original treatment. In general, we follow [Wolhre et al.| and 
set <))pah = 0.5. 

One obvious question that arises at this point is whether 
it is valid to take a prescription for the photoelectric heat¬ 
ing rate that was designed primarily to model the behavior 
of the diffuse ISM and to apply it to a dense PDR. After 
all, the interstellar radiation held seen by a representative 
patch of the diffuse ISM has a spectral shape which differs 
signihcantly from that of a blackbody. Therefore, before we 


can apply the Bakes & Tielens (1994 1 prescription for photo¬ 


electric heating, we have to make sure that the assumptions 
that they made in the derivation of Ppe are not too strongly 
violated for our radiative transfer method. There are two 
main concerns: the spectral shape of the ISRF is not that of 
a blackbody, and absorption by dust is not the only attenu¬ 
ating process in the relevant energy range. 

We quantify the deviation from the ISRF by taking a 
reference wavelength of 1110 A (i.e. 11.17 eV) and normal¬ 
ize blackbody spectra at different effective temperatures so 
that they all have the same energy density at this refer¬ 
ence wavelength (see Figure]^. For effective temperatures 
Tefi = 20000-50000 K that are relevant for massive stars, 
we find that the energy densities of our normalized spectra 
at the wavelengths of interest typically differ by no more 


than a factor of three from the Draine (1978 1 model for the 


ISRF. Bakes & Tielens (19941 calculate the heating rate for 
a radiation field with the shape of a Teff = 3 x 10"^ K black¬ 
body and compare this to the results they obtain for the 
Draine ISRF, normalized so that Go is the same for both 
spectra. They find that in this case, the difference in spec¬ 
tral shape makes around a 25% difference to the photoelec¬ 
tric heating rate. A spectrum with effective temperature of 
Tefi = 2 X 10'* K nearly reproduces the Draine field and so 
in this case, we would expect the error to be much smaller, 
around a few percent. For higher effective temperatures of 
4 X 10* K and 5 X 10* K, the shape of the spectrum does 
not change drastically from the reference 3 x 10* K black¬ 
body, and so although the error will be larger, we would still 
expect it to be less than 50%. We are therefore justified in 
using the same efficiency parameter e as in [Wolfire et al.| 
(20031, as any errors introduced by the difference in spec¬ 


tra shape are relatively small and are probably dwarfed by 
the uncertainties stemming from other features of the dust 


physics, such as the PAH abundance in dense clouds or the 
exact geometric shape of the dust grains. 

As a star evolves on the main sequence and beyond, 
its spectrum changes. We use non-rotating, non-magnetized 
solar metallicity tracks for typical massive stars taken from 
Ekstrom et al. (|2012[) and [Georgy et ah (20121 to illustrate 


that these stars stay inside the explored range of Tefi (see the 
right hand side of Figure]^. A massive star enters the Wolf- 
Rayet phase at the end of its lifetime t*, at which point its 
effective temperature varies rapidly. For this stage, our mod¬ 
eling for the photoelectric heating efficiency breaks down as 


equation (491 is no longer applicable. Fortunately, only a 


short period of about 10 to 15% of t* is spent in this phase. 

The other issue that we need to address is how to prop¬ 
erly account for the fact that photons in the energy bin 
can be absorbed by either dust or H 2 . In principle, we should 
treat the photons in this energy bin in the same careful way 
that we treat the photons in the E \^,2 energy bin, by first cal¬ 
culating the total number that are absorbed and then parti¬ 
tioning them between H 2 and dust. However, in practice this 
proves to be unnecessary. The reason is that in the regime 
where both H 2 photodissociation and photoelectric heat¬ 
ing are important, the effective H 2 photodissociation cross- 
section (Th 2 is orders of magnitude larger than the effective 
dust absorption cross-section, which for typical ISM dust 
peaks at around 700 A with a value of Cd « 3 x 10 “^* cm^/H 
and then drops off rapidly at both larger and smaller wave¬ 
lengths ( Draine|2003 Weingartner k, Draine|[200T l. There¬ 
fore, in these conditions H 2 absorption dominates. This 
means that rather than working out the number of photons 
absorbed by the combined effects of H 2 absorption and dust 
absorption and then partitioning these photons up between 
the two absorbers, we can instead safely split up the calcu¬ 
lation into two steps: we first compute how many photons 
are absorbed by H 2 and remove these from the i?ii ,2 energy 
bin, and only then use the remaining flux to estimate Go 
and Fpe. 

In practice, it is convenient to use the incident energy 
flux, Tpe in erg s“*, in place of the photon flux when we es¬ 
timate the local stellar intensity Go in units of the Habing 
field. This consists of two contributions. The first term treats 
photons with energies in the range of 5.6 eV to 11.2 eV. 
These photons are unable to electronically excite H 2 from 
the ground state and so are only absorbed by dust. The 
second term adds the contribution from photons in the en¬ 
ergy range 11.2 eV to 13.6 eV that have not already been 
absorbed by molecular hydrogen. We therefore have: 

Epe = T5.6(r = 0)e-^^v + 

(51) 

(/shdA''ll. 2 e ’"'‘/At — (1 + /pump) fcuv)(7^11.2) 

where = 0 ) is the energy flux at the stellar surface in 

the energy range from 5.6 eV to 11.2 eV, which we evenly 
distribute over the initial rays, and 7 = 2.5, taken from 


Bergin et al. (20041. To determine the FUV energy flux, we 


first determine the FUV photon flux entering the cell and 
then subtract from this the photon flux absorbed by H 2 in 
the cell. Finally, to convert the result from a photon flux to 
an energy flux, we multiply by the average energy (i 5 ii. 2 ) of 
the photons in the E \\,2 energy bin. 

The energy flux = 0) is carried along the ray as a 

ray property and attenuation is always calculated from this 
value. The process is not photon-conservative as we do not 
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Figure 3. Left panel: blackbody spectra for different effective temperatures as well as the Habing and Draine ISRFs. All spectra are 
normalized to their energy density at 1110 A. The black dashed line is the Habing ISRF and the black solid line is the Draine ISRF. 
Colored solid lines correspond to effective temperatures of 2 X 10^ (blue), 3 X 10^ (green), 4 X 10^ (orange) and 5 X 10"^ (red) Kelvin. 
Right panel, the evolution of the effective temperatures over the lifetime t* of typical massive stars with tracks taken from |Ekstrom et al.| 
l |2012[ | and [Georgy et aL] ( |201^ . Here the colors denote the mass of the star, where the more massive stars have a higher blackbody 
temperature. The masses are, from top to bottom, 60, 25, 15 and 8 Mq. 


equate the extinction of the photo-electric heating photons 
with the eventual number of photons that heat the surround¬ 
ing gas, but this is unlikely to introduce a signihcant error 
into our calculation of the heating rate. 

After the ray-tracing step is complete and the energy 
flux of all contributing rays in each cell is added up, Ape has 
to be expressed in units of the integrated Habing field as 
input for equation ( |49| |. We convert Ape to an approximate 
flux per unit area by dividing it by the area of one face of the 
cell, cP. We then divide the resulting energy flux per unit 
area by the fiducial value of Gconv = 1.6 x 10“® erg s“^ cm“^ 


|Habing|1968 1 to obtain Go, 


Go = 


d^Gc 


(52) 


Note that this conversion to Go naturally accounts for adap¬ 
tively refined cells with different physical sizes d. 


2.2.6 Solving the energy equation 

In order to solve for the evolution of the internal energy, we 
need to know the heating rates due to photoionization, H 2 
photodissociation, the UV pumping of H 2 and photoelectric 
heating, which are computed as described in the previous 
sections. In addition, we need to know the radiative cooling 
rate of the gas, and must also account for any other chemi¬ 
cal heating or cooling terms (e.g. H 2 formation heating, H"*" 
recombination cooling). To compute these additional terms, 
we use a modified version of the heating and cooling function 
introduced in jGlover et al.| ( |2010| | and updated in [Glover fcj 
Clark (2012a I. This accounts for all of the important chem¬ 


ical heating and cooling terms, as well as radiative cooling 
from the fine structure lines of C*" and O, the rotational 
and vibrational lines of H 2 and CO, the electronic transi¬ 


tions of atomic hydrogen (“Lyman-a” cooling) and thermal 
emission from dust. We have supplemented this by includ¬ 
ing the additional forbidden and semi-forbidden transitions 
of C^, O, O'*" and N"*" that are summarized in Table 9 of 


Hollenbach & McKee (19891, as these become important 
coolants close to T ~ 10^ K and hence make an appre¬ 
ciable difference to the equilibrium temperature of our Hll 
regions. As we do not track O'*' and N"*" explicitly in our 
chemical network, we make the simplifying assumption that 
* 0 + /®o = a;N+ /®N = ®H+ /s^H. Finally, we also include cool¬ 
ing from the electronic transitions of atomic helium and from 


metals, using the rates tabulated by Gnat & Ferland (20121. 


These rates assume that the gas is in collisional ionization 
equilibrium, which is not the case within our Hll regions, 
since the latter are dominated by photoionization. However, 
these processes become important in comparison to Lyman- 
a cooling only aX T 10"^ K, and hence do not strongly 
affect the thermal evolution of the gas in our simulations. 
Nevertheless, these high temperature cooling processes are 
needed for a complete treatment of the ISM temperature 
structure and will be neccessary for future extensions of the 
simulation code. 

Putting all of these pieces together, we can write the 
final net heating/cooling rate as 


Ftot = (F — A)g10 — AhM89 — Agfi2 

+ Fuv + Fpe -|- Fion + Fdis, 


(53) 


where (F — A)gio denotes the heating and cooling function 
from Glover et al. (20101, Ahm89 corresponds to the forbid¬ 


den and semi-for bidden line cooling from [Hollenbach &: Mc-[ 
|Kee| ( [T98^ and Agfi 2 corresponds to the high-temperature 
helium and metal-line cooling from Gnat & Ferland 120121. 
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2.2.7 Missing physics 

In our radiative transfer scheme we do not treat any scat¬ 
tering processes or diffuse radiation. This leads to shadows, 
cast by optically thick gas, that are too sharp. 

In addition, radiation pressure in the form of momentum im¬ 
parted on H nuclei is also not taken into account. It could be 


included in a similar fashion as in Wise & Abel (2011 1 , who 


modeled dust-free primordial gas, but as we aim for a dusty 
present-day ISM, an isolated treatment acting only on hy¬ 
drogen is incomplete. Ideally, an effective radiation pressure 
based on the complete spectrum, dust composition and size 
distribution should be implemented. Furthermore, in harsh 
environments such as in close proximity to massive stars the 
assumption that dust is tightly coupled to gas would have 
to be reevaluated. 

We do not include any ionization of other elements be¬ 
sides hydrogen. Energetically, only helium has a small effect 
on the overall energy deposition as it makes up about 10 % 
of the ISM in abundance. Other elements, such as oxygen 
or nitrogen, only exist in trace amounts unable to capture 
any significant amount of radiation. We can roughly esti¬ 
mate the upper limit of the thermal energy deposition from 
ionization of helium for the most extreme blackbody tem¬ 
perature of Teff = 50000 K we consider here, by considering 
an ISM that consists only of helium and one that is pure 
hydrogen. Integration of the blackbody spectrum with an 


He ionization cross-section taken from Verner et al. (1996 1 


yields a total possible energy deposition ten times less in the 
case of the helium ISM in comparison to the hydrogen one. 
Taking into account that only around one-tenth of the inter¬ 
stellar medium is made up of helium, this roughly amounts 
to a percent effect in additional energy converted from the 
radiation field. 


3 TESTS 

In this section we study the influence of spatial and tem¬ 
poral resolution on the expansion of ionization and photo¬ 
dissociation fronts. Our testing strategy follows the standard 
approach of using simplified setups to check the included 
physics piece by piece and then in combination. The first 
three tests only include atomic hydrogen ionization, in a 
static and dynamic density field, with known analytical so¬ 
lutions for the expansion of the ionization front in its initial 
(R-type) and later stage (D-type). 

The fourth test checks the non-ionizing radiation cou¬ 
pling to the chemical network in a test designed for pho¬ 
todissociation region (PDR) codes. It first probes the chem¬ 
ical state of the gas at fixed temperature, then the thermal 
state of the gas is allowed to change. The last two test cases 
are of a qualitative nature, where we include all modeled 
physics. Here, we check for any apparent numerical artifacts 
or unexpected behavior in the radiation-gas coupling. 


3.1 R-type ionization front expansion 

The ionization front (I-front) expansion velocity during its 
R-type phase is much greater than the sound speed of the 
ambient neutral gas. Hence, its supersonic expansion leaves 


no time for the photo-heated, over-pressured gas to com¬ 
press the surrounding medium. Numerically, simulating the 
R-type phase tests the ray-tracing and rate calculation al¬ 
gorithm as well as the chemical network, although for sim¬ 
plicity we ignore H 2 in this test and simply consider the 
ionization of atomic hydrogen. 

|Stromgren| ( |1939[ ) derived the size of an ionized (Hll) 
region by considering the equilibrium between photoioniza¬ 
tion and recombination. In a uniform density gas of pure hy¬ 
drogen, the ionized region is spherical (a Stromgren sphere) 
with radius 


R. 


f i Q y 

\47r n^QB ) 


(54) 


where n is the hydrogen nuclei number density of the am¬ 
bient medium, Q — dWon/dt the rate at which ionizing 
photons are produced by the source, and qb is the case B 
hydrogen recombination rate, where recombination to the 
ground state is excluded, because the emitted photon is as¬ 
sumed to be quickly absorbed again by a hydrogen atom in 
the vicinity with no net change in the ionization fraction. 
The time evolution of the radius of the ionization front as 
it approaches this equilibrium solution can be recovered by 
considering the possible fates for the ionizing photons. In a 
dust-free gas of pure hydrogen, all of these photons are either 
absorbed within the Hll region, where the ionization state 
has to be maintained against continuous recombination, or 
propagate through it and reach the shell-like ionization front 
at position r-,, where they ionize the gas, enlarging the Hll 
region. If dWon photons reach the I-front, then they ionize 
a thin shell with thickness dri, where dNion ~ dnnrfdri. It 
therefore follows that 


diVi, 

dt 


2 dr, 4 32 

= 47 rri n— - 1 —yrr; n as 

dt 3 


(55) 


where the second term on the right-hand side accounts for 
recombinations within the Hll region, which we have as¬ 
sumed is almost completely ionized. Rearranging this ex¬ 
pression and solving for ri(t) then yields 

n(t) =Re(l-e-‘/‘"‘=)^/^ (56) 


where tree = (obu)”^ is the recombination time. 

For our test, we set up a simulation in which we fix 
qb to the constant value qb = 2.59 x 10“^® cm^ s“^ and 
consider gas that has a quasi-isothermal equation of state, 
with 7 = 1.0001. The atomic hydrogen number density is 
set to n = 10 “® cm“®, and we consider a simulation do¬ 
main that extends from 0 to 6.4 kpc in the x-, y- and 2 ; 
directions. The source emits Q = 10'*® photons per second, 
and we assume that these photons are monochromatic, with 
energy E = 13.6 eV, the ionization threshold of hydrogen. 
This choice means that there is essentially no photo-heating, 
and also allows us to consider a fixed photoionization cross 
section, cth = 6.3 x 10“*® cm®. We capture an octant of the 
spherical I-front expansion by positioning the source at the 
origin. For the purposes of this simple test, we neglect the 
effects of dust and consider no chemical processes other than 
photoionization and radiative recombination. 

Figures |4||^ show the impact of temporal and spatial 
resolution on the I-front expansion. In these figures, we de¬ 
fine the position of the ionization front as the point where 
the neutral hydrogen fraction drops below 50%. 
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Figure 5. Impact of spatial resolution on the expansion of the 
I-front during its R-type phase. The thin solid line shows the 
analytical solution, while the four thick lines show simulations 
with 64^ (dotted), 128® (dash-dotted), 256® (dashed) and 512® 
(solid) cells. In all four simulations, we set /jj = 0.1. 


Figure 4. Top: R-type ionization front expansion in a uniform 
medium. The solid line shows the analytical solution, the dotted 
line is without a timestep limiter, the dash-dotted line shows the 
expansion with the change per timestep in the atomic hydrogen 
fraction limited to 10% (i.e. /h = 0 .1), and the dashed line shows 
the case where the change is limited to 1%. Bottom: Relative error 
in our solution compared to the analytic solution. The line styles 
are the same as in the top panel, but we omit the case without 
the timestep limiter. 


Figure [^demonstrates the necessity for a timestep lim¬ 
iter, as without it, we systematically underestimate the I- 
front radius, in this case by as much as 50% at t ~ 0.3 tree. 
This occurs because in our calculation of dr for each grid 
cell, we implicitly assume that the hydrogen number den¬ 
sity does not vary significantly during the timestep. If we 
do not employ a limiter, then this assumption is not al¬ 
ways justihed and can lead us to overestimate the number 
of photons that are absorbed in each cell. However, we also 
see that if we limit the change of neutral hydrogen to 10 % 
with /h = 0 . 1 , the expansion matches the analytical calcu¬ 
lation to within an error of only a few percent. Decreasing 
/h further, to /h = 0 . 01 , improves the solution even more, 
but the improvement is relatively small in view of the large 
computational cost. 

The difference between our numerical solution and the 
analytical solution is largest at early times, since at this 
point the neutral hydrogen fraction is changing rapidly in 
the intense radiation field close to the source. Our timestep 
limiter can detect this and reduce the size of the timestep 
to compensate, but only after we have already taken one 
step with a too-large timestep. It therefore reduces the error 
only on subsequent timesteps, which for large /h is not quite 
enough to recover the proper I-front expansion. 

As far as sensitivity to the spatial resolution is con¬ 
cerned, we see from Figures [^ and [^ that if we set /h = 0.1, 
we detect virtually no change in the expansion behavior as 
we increase the spatial resolution from 64® to 512®. The 


main effect of the change in the resolution is to alter the 
thickness of the I-front, since it has a minimum width of 
three resolution elements. In addition, the initial expansion 
is slightly influenced by the fact that at high resolution, the 
cells have smaller volumes and hence will be ionized more 
quickly compared to cells in more coarsely resolved simula¬ 
tions. This means that high resolution simulations are more 
sensitive to our choice of /h than low resolution simulations, 
since the initial change in nu in the central cells occurs more 
rapidly as we decrease the size of the cells. Note, however, 
that the error introduced by this only really affects the evo¬ 
lution of the I-front during the first few simulation steps. 

Finally, Figure also demonstrates that our ray-tracing 
scheme properly captures the spherical nature of the I-front 
without introducing any artifacts due to the underlying 
Cartesian mesh. 


3.2 Ionization front expansion in an r ® density 
profile 


Hll regions in simulations of star formation usually do not 
expand into a uniform medium at rest. Instead, the ioniza¬ 
tion front usually expands along a density gradient from a 
dense core. We test the applicability of our radiative trans¬ 
fer scheme in this case by assuming a density profile of the 
form 


n{r) 


ric if r ^ rc 

nc{r/rc)~‘^ if r > r^, 


(57) 


with radius r from the core center, and where ric is the 
hydrogen number density in the core and is the size of 
the constant density center of the core. In such a density 
profile, equation (551 becomes 


~ 47i'rfn(ri)‘^ -f 47rQB f n{r)^C dr. (58) 

dt dt Jq 

We are interested in the case where the I-front has already 
left the central dense clump, i.e. ri > r^. If we fix the ex¬ 
ponent of the density profile to be cu = 2 , then it is easy to 











14 


C. Baczynski, S. C. O. Glover, R. S. Klessen 


64 ^ 128 ^ 256 ^ 512 ^ ®hii 



X (kpc) X (kpc) X (kpc) X (kpc) 


Figure 6 . Slices through the simulation domain of the R-type test for runs with 64, 128, 256 and 512 cells per side taken at f 130 Myr. 


show that 


dri Q 4 , UcrluB 

Id ~ RRRkI ~ . 


(59) 


where Q = dNion/dt and Vi = dri/df. If we now choose Q 
so that the first two terms on the right hand side of this 
eqnation snm to zero, then this eqnation has the simple an¬ 
alytical solution (Franco et al. 1990 [Mellema et al. 2006 
[Whalen fc Norman|2006P 


n(t) = rRl + 2 trie 


(60) 


where we have taken t = 0 to be the time at which the I-front 
escapes from the constant density portion of the profile, i.e. 
when n = Tc. 

To test whether our code can reproduce this analytical 
solution, we make the same approximations as before. We 
set 7 = 1.0001, consider monochromatic photons with E = 
13.6 eV, so that ch = 6.3 x 10“'^® cm^ and there is no 
photo-heating, and fix ob to the constant value qb = 2.59 x 
10“'^® cm® s“^. We set Q = 10"^® photons s“® as before, and 
consider a central density ric = 100cm“®. With these values, 
we then need to set — 1.99 pc in order to ensure that the 


first two terms in equation (591 vanish. 


We consider a simulation domain that extends from 0- 
6.5 pc in each dimension, and center our spherical density 
profile within this domain. The source is placed at the center 
of the density profile. As we want to test the behavior of the 
I-front in the regime where > r^, we assume that all of 
the gas at r ^ rc is already ionized. 

The main issue tested with this setup is whether our 
timestep limiter is able to cope with the sustained rapid 
change of the neutral hydrogen fraction during the expan¬ 
sion of an I-front down a steep density gradient. The top 
panel of Figure shows the ratio of the ionization rate to 
the number density of H atoms at the I-front for both the 
uniform density and the power-law profile cases. It is clear 
that this quantity varies much more rapidly with the power- 
law density profile n = const., particularly when t < tree. 

The middle panel in Figure]^ shows the position of the 
I-front as a function of time in simulations with a spatial res¬ 
olution of 256® grid cells and with /h = 0.1 (dotted line) and 
/h = 0.01 (dashed line). The analytical solution is shown as 




Vhec 

Figure 7. Top\ Ionizations per second normalized to the atomic 
hydrogen number density at the I-front position. The dashed 
curve corresponds to the case of uniform density (tih = const.), 
while the solid line corresponds to the density profile given by 
equation l |57[ l. Middle: comparison of the numerical to the ana¬ 
lytical solution for the I-front expansion with different timestep 
limiters: The dotted line shows the expansion with fn limited to 
0.1, while the dashed line shows the result for /h = 0.01. Bot¬ 
tom: The fractional deviation of the numerical result from the 
analytical solution. 
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the solid line. We find that limiting the change in the atomic 
hydrogen fraction in each cell to 10 % per timestep allows 
us to reproduce the analytical solution to within around 5% 
(see the bottom panel of Figure]^. The error becomes larger 
as the medium rarefies and the change in ionization fraction 
in a single cell increases in each timestep. Using a stricter 
limit with /h = 0.01 allows us to reduce the error to less than 
a percent, but significantly increases the computational cost. 


3.3 D-type ionization front expansion 


Once the expansion speed of the ionization front drops below 
the speed of sound in the hot, over-pressured ionized gas, 
the shock-front traveling on top of the I-front detaches and 
propagates supersonically into the cold surrounding ambient 
medium. The shock starts to sweep up a dense shell, and the 
I-front enters its D-type phase. If the neutral gas is very cold, 
then we can neglect its thermal pressure. In this case, the 
rate of change in the momentum of the dense shell located 
ahead of the I-front is equal to the total force exerted on the 
shell by the pressure of the hot gas, i.e. 


dt 


{rriiVi) = 4:nrfpi, 


(61) 


where pi is the pressure of the ionized gas. Note that we have 
also assumed here that gravity is unimportant and that the 
radiation pressure of the photons can be ignored. 

Let us now consider an ionization front expanding into 
a constant density medium with p — po- If we assume that 
the gas that was initially in the region occupied by the Hll 
region has all been swept up into the shell, then the shell 
mass is simply 

mi = —ripo. (62) 

If the shell is thin, then the pressure pi of the ionized gas 
is balanced by the ram pressure exerted by the neutral 
medium, pi = picf = pov^, where Vs is the shock expan¬ 
sion velocity, a is the sound speed in the ionized gas and pi 
the density of the ionized gas. 

If we also assume that the Hll region is in ionization 
equilibrium, so that Q = ( 47 r/ 3 )ni oerf, i.e. all of the ion¬ 
izing photons are used up counteracting the effects of re¬ 
combination, then we can express the density in the ionized 
region as 


Pi = Po 



(63) 


where Ra is the Stromgren radius. With these assumptions, 
equation (611 reads 


d/'4 3 .\ 2 2 

dt V 3^’'“ ^ 


3/2 


(64) 


Hosokawa & Inutsuka (2006 1 solve this equation and arrive 


at the expansion law for a D-type ionization front 

4/7 

n«) = a(i + yig' 


(65) 


(Note that the more commonly known solution by Spitzer 
(1978 1 omits the factor of ^4/3 from the second term inside 


the parentheses). Eventually, the expansion of the I-front 


should come to rest once the pressure exerted by the ionized 
region is balanced by the pressure of the ambient medium 
( Raga et al.|2012 |. 

We use the D-type phase of the ionization front ex¬ 
pansion as a test for the hydrodynamical response to the 
over-pressured Hll region and the continuous driving of the 
shell by ionizing radiation. The simulation setup consists of 
a source at the origin that emits Q = 10^® ionizing photons 
per second, with 2 eV per photon in deposited heating en¬ 
ergy. The hydrogen ionization cross-section is again fixed at 
the threshold value of cth = 6.3 x 10“^® cm^. 

We set the atomic hydrogen number density to nn = 
1000 cm“®, and also set 7 = 1.66667 and as = 2.59 x 
10“^^ cm® s“®. The ambient temperature of the neutral gas 
is set to 100 K. The extent of the simulation domain is cho¬ 
sen to be ±64 pc. To balance the photoionization heating, 
we consider only Lyman-a cooling; the rates of the other 
cooling and heating terms in our model are set to zero. We 
perform simulations with spatial resolutions of 64®, 128®, 
256® and 512®. 

With this problem setup, the initial Stromgren radius is 
Rs ~ 0.68 pc. This is resolved with multiple grid cells only 
in our 512® simulation. Figure shows the expansion be¬ 
havior of the I-front for different spatial resolutions. When 
the initial Stromgren sphere is unresolved, the shock does 
not form until a fully ionized and over-pressured region has 
been created. This is a purely numerical effect and is due to 
the fact that for low resolution, the ionizing photons are di¬ 
luted over a large volume, only slowly heating and changing 
the ionization fraction. Once a minimally sized Hll region of 
one or two cells is created, the D-type expansion progresses 
as described by the analytical expansion law. Only if the 
Stromgren sphere is at least marginally resolved do we see 
no delay or a negligible delay (see Figure]^. We correct for 
this time offset by replacing the evolution time t in equation 
( |65[ ) with t' = t — toff, where toff is the time it takes to ionize 
the minimum volume of one cell. 

We use the speed of sound calculated from the average 
temperature in the ionized gas Ci = yfesTavg//!, with the 
Boltzmann constant ks and mean molecular weight p. At 
later times deviations from the analytical solution emerge 
that are a consequence of the fact that the ionized gas is 
not isothermal. We find that in practice, the lowest temper¬ 
ature in the ionized gas, is located close to the I-front 
(see Figure]^, and that at late times, the expansion of the 
ionization from follows the analytical solution that we ob¬ 
tain if we set T — Tmin rather than T — Tavg. In general, 
the difference between our numerical solution and the an¬ 
alytical solution is around 5%, apart from the initial rapid 
expansion phase, where for an unresolved Stromgren sphere 
the error reaches up to 30%. 

Figure 1^ shows the impact of the spatial resolution on 
the shape of the shell, as well as the temperature and density 
structure of the Hll region. We see that the temperature 
varies by around a factor of two as we move from the central 
region around the star to the edge of the dense shell. This 
temperature variation is a consequence of spatial variations 
in the heating rate and the gas density (which directly affects 
the cooling rate). We note, however, that the effect is large 
here because we are not including the contribution of metals 
to the total cooling rate. If we include the additional metal¬ 
line cooling, we find instead that the temperature of the 
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Figure 8. Expansion of the D-type ionization front in our test setup with different spatial resolutions. The dotted, dashed and dash- 
dotted lines show analytical solutions with the sound speed c, calculated from the maximum (Tmax), average (Tavg) and minimum (Tmin) 
temperature of gas with an ionized hydrogen fraction of ^ 0.99. The solid lines plot the expansion obtained from the numerical 
simulation after correcting for the offset time needed to generate the initial ionized and over-pressured Hll region. The bottom panels 
show the relative errors calculated from the analytical solutions where we use Tavg (dashed) and Tmi^ (dash-dotted). 


ionized gas evens out at around 8-9 x 10^ K with a negligible 
temperature gradient. 

Comparing the ionized fraction with the density slice re¬ 
veals that the shock has a direction-dependent expansion ve¬ 
locity: it propagates slightly faster in directions aligned with 
the Cartesian mesh and slightly slower in the other direc¬ 
tions. This is a well known effect for spherical shock waves on 
Cartesian grids and can largely be eliminated by increasing 
the resolution of the simulation, as Figure [^demonstrates. 


3.4 Photon dominated region test 


We follow the test setup outlined in the code comparison 
project in [Rollig et al.| ( |2007| ). The test was designed to 
compare dedicated PDR codes that obtain the chemical and 
thermal structure of a molecular cloud illuminated by the 
interstellar radiation field. The final, equilibrated state is 
calculated and compared, where any time evolution is omit¬ 
ted. The chemical network used in their test case consists of 
31 different chemical species and hence is far more extensive 
than the simplified network used in Fervent. Despite this, 
it is interesting to compare how well our code can reproduce 
the chemical structure of the two main molecular species, 
H 2 and CO, and the thermal structure of the PDR. 

Table[^shows the elemental composition of the molecu¬ 
lar cloud, as well as the test parameters. In our simulations, 
the chemical abundance of atomic oxygen is tracked using a 
conservation law: we know the total elemental abundance of 
oxygen relative to hydrogen (Ao) and also how much oxy¬ 
gen is incorporated into CO. The fractional abundance of 
atomic oxygen therefore follows as xo = Aq — xco- Our 
chemical model does not include neutral carbon, and so we 
assume that any carbon not incorporated into CO is present 
as C"*". 


The tests in [Rollig et al.| ([2007| are carried out using a 
UV radiation field that has the spectral shape of the |Draine[ 
119781 field, and a strength parameterized by x, where X = 1 
corresponds to the Draine field. This UV radiation field pro¬ 
duces an H 2 photodissociation rate in unshielded gas that is 


Table 2. Parameters of the PDR test 


Symbol 

Value 

Quantity 

d-He 

0.1 

He abundance 

^0 

3 X 10-"^ 

0 abundance 

do 

1 X 10-^ 

C abundance 

kcT 

5 X 10“^'^ s-i 

cosmic ray ionization rate 

dv 

6.289 X 10-227Vh, total 

visual extinction 

'^d 

3.02Av 

dust attenuation 

fcuv 

5 X 10-“(x/10) s-i 

H2 dissociation rate 

i?H2 

3 X 10-i®tV 2 cm^s-i 

H2 formation rate 

Tgas 

50 K 

isothermal gas temperature 

Tdust 

20 K 

isothermal dust temperature 

■fd/g 

1.0 

dust to gas ratio 


fcuv = 5 X 10“^''(x/10) s“^, as listed in the Table. However, 
in Fervent we do not directly specify fcuv, and so we must 
instead calculate the photon flux in the 11.2-13.6 eV energy 
bin that we require in order to produce the same dissocia¬ 
tion rate. We can write the energy density of the adopted 
UV radiation field as | Draine|1978 Draine & Bertoldi 19961 


where A is the wavelength in A, A 3 = A /1000 A, and ux is 
the energy density per unit wavelength, measured in units of 
erg cm“^ A“^. Three quantities are needed as input for our 
scheme: the incident photon flux Afii. 2 (r = 0 ) in the E 11.2 
energy bin 


A /■mo Xu 

NuRr = 0) = - / -^dA = 1.06 x lo" xA, (67) 

where A is the area of the illuminated face of our simulation 
domain and the factor of one half enters because we are con¬ 
sidering a semi-infinite plane parallel slab which is exposed 
to radiation from a total solid angle of only 27r steradians; 
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Figure 9. Slices of the D-type expansion test through the origin showing different hydrodynamical quantities after tree ~ 81000 
recombination times or 9.5 Myr. Spatial resolutions of 2 pc (64^), 1 pc (128^), 0.5 pc (256®) and 0.25 pc (5123) are shown. 


the average energy per UV photon in erg, 


rlllO rl 

{E 11 . 2 } = / cuxdx/ j 
J912 ' J91 


Ama 


dA 


/912 
= 1.93 X 10“ 


( 68 ) 


where we have to express the speed of light in A s”^; and 
finally, the total energy flux in the i? 5.6 energy bin at the 
cloud surface, measured in units of erg s“^cm“3 

1 c2200 

Fi.eir = 0) = - / cuxdX = 1.22 x 10"^ x- (69) 
^ iiiio 

This last quantity is needed in order for us to be able to 
account for the effects of photoelectric heating when deter¬ 
mining the thermal structure of the PDR. 

The simulation domain is an elongated box with a ratio 
of (x, y) : z oi 1 ■. 10, with the physical extent set to give 
a range in visual extinction depending on the gas density 
used in each test. The ray-tracing is simplihed to cast rays 
without the HEALPix based formalism. Instead, the initial 
rays are created at the centers of the faces at the minimum 


X boundary, which we assume to be illuminated by the ISRF. 
They propagate in the positive x direction, where they are 
split into four child rays if they encounter a more highly 
resolved region, one for each child cell (see Appendix [A] for 
a brief description of the mesh structure). 


The PDR comparison project described in|Rollig et al.| 


(20071 consists of eight tests. Four of these are isothermal 
with fixed gas and dust temperatures (denoted FI to F4). 
The other four are carried out with the same density and 
ISRF strengths as the isothermal tests, but the gas and dust 
temperatures are allowed to evolve self-consistently with the 
chemical state. These non-isothermal tests are denoted as 
V1-V4. The tests explore two different densities and two 
different radiation held strengths: a cloud with a hydrogen 
number density of nu = 10^ cm~3 illuminated by an ISRF 
with X = 10 or X = 10® (tests Fl-2, Vl-2) and a cloud with a 
hydrogen number density of nu ~ 10®'® cm”® with the same 
variation in the ISRF held strengths (tests F3-4, V3-4). 


In the isothermal test cases, we use an adiabatic index 
of 7 = 1 . 0001 , while in the non-isothermal tests we adopt the 


./§] (cy)So, 
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standard 7 = 1.6667. The extent of the simulation domain is 
chosen so as to provide a sufficiently high spatial resolution 
to sample small values of visual extinction Av- For tests Fl- 
2 and Vl-2, the smallest cell size is set to 2.44 x 10^'* cm 
and the simulation domain extends to ±2 x 10 ^® cm in the 
X and y directions and 4 x 10^® cm in the 2 : direction. For 
tests F3 and V3, the bounding box size is ±0.5 x 10^® cm 
in X and y and 10 ®"^ cm in the z direction, meaning that the 
most resolved cell has a size of 6.1 x 10® cm. For tests F4 
and V4 the box extends to ±0.5 x 10®® cm in x and y and 
10 ®®” cm in the 2 direction, with the smallest cell having a 
size of 2.44 x 10®® cm. 

In contrast to many dedicated PDR codes, our radiative 
transfer scheme has to be evolved explicitly in time until 
the chemical network reaches its equilibrium. During this 
evolution, we do not follow the hydrodynamical response 
of the gas to any changes in the thermal pressure, i.e. the 
density field remains fixed in space and constant. 

Figure shows the results that we obtain for the four 
isothermal tests. For comparison, we also plot the results 
from the eight PDR codes compared in Rollig et al. (20071, 
focusing on the surface densities of H, H 2 , C"’" and CO. We 
see from the Figure that Fervent does a good job of repro¬ 
ducing the transition from atomic hydrogen to molecular hy¬ 
drogen at the edge of the cloud. The PDR codes compared 
in Rollig et al. (20071 disagree somewhat on the details of 


this transition, and our results lie well within this spread. 

On the other hand, it is clear from the lower panels 
of Figure that we do not reproduce the behavior of the 
CO or surface densities. The main reason for this is 
our assumption that we can obtain the CO dissociation rate 
from the H 2 dissociation rate simply by applying a conver¬ 
sion factor. In reality, this is accurate only in optically thin 
gas. Once the gas starts to become optically thick, H 2 self¬ 
shielding plays a major role in decreasing the H 2 photodis¬ 
sociation rate (Draine & Bertoldi 19961. However, the anal¬ 
ogous process for CO is much less effective, on the grounds 
of the low abundance of CO relative to H 2 , and in any event 
does not have the same functional dependence on column 
density as the H 2 self-shielding correction (see e.g.|Lee et al.| 


1996). By directly coupling the CO and H 2 photodissocia¬ 


tion rates, we therefore overestimate the rate at which the 
CO photodissociation rate falls off as we move into the cloud, 
and hence overestimate the CO surface density. In addition, 
our extremely simplified treatment of CO formation is also 
known to overestimate the CO formation rate relative to 
more accurate models ( Glover fc Clark|2012a I. We therefore 
see that we cannot use Fervent in its current form for pre¬ 
dicting accurate CO abundances. However, we stress that 
our goal when adding the current simplified treatment of 
CO chemistry was simply to be able to approximately dis¬ 
tinguish between regions of the ISM that are cold and CO- 
bright and ones which are warm and CO-faint when per¬ 
forming large-scale simulations (see e.g. Walch et al.|2014 |, 
and for this particular purpose, our current treatment is ad¬ 
equate. In the future, we intend to improve on Fervent’s 
treatment of CO by adding an additional energy bin specif¬ 
ically to treat CO photodissociation, and by using a more 
sophisticated treatment of CO formation, such as the model 
of Nelson & Langer (1999), but this is outside the scope of 
the present work. 

Finally, we also see from Figure that Fervent over¬ 


estimates the amount of ionized carbon in the PDR by 
roughly an order of magnitude. This is another consequence 
of our reduced chemical model: since we do not include neu¬ 
tral atomic carbon, we mis-classify as C^ some carbon that 
should in fact be present as C. Again, we plan to address 
this in the future work by improving the model of the carbon 
chemistry. 

In Figure EH we show the results from the non- 
isothermal tests. We again find a good agreement in most 
cases between our results for the H 2 and H surface densities 
and those computed by the PDR codes compared in |Rolli^ 


et al. (2007), particularly for the lower density PDRs. The 


gas and dust temperature structure itself lies within the very 
large spread obtained from the dedicated codes, except for 
test V2 where we find a temperature that is around a factor 
of two lower than in Rollig et al. (2007). We note that one 
factor that may contribute to this is the fact that in our 
treatment of fine structure cooling from atomic oxygen, we 
use collisional excitation rate coefficients for O-H collisions 
taken from Abrahamsson et al. (2007). These are a factor of 
two to three higher than the older Launay & RouefI (1977) 
values available at the time that the [Rollig et al.| code com¬ 
parison was carried out, and so in conditions where atomic 
oxygen cooling dominates (warm, dense, atomic gas), we 
would naturally expect to recover lower temperatures than 
in these calculations. Finally, we again see that we do not do 
a good job of reproducing the behavior of the CO and C"*" 
surface densities, for the same reasons as explained above. 


3.5 Isolated source in a dense molecular medium 

In this test, we again examine the expansion of an Hll re¬ 
gion around a point source embedded in a uniform density 
gas. The difference is that in this case, we take the gas to 
be fully molecular initially, and we account for the PDR 
that forms ahead of the I-front. This is a complex multi¬ 
physics problem, where numerical simulations are actually 
the best tool for exploring the multiple combined heating 
and dissociation processes. No simple anal 3 ±ical solutions or 
published test results are available for comparison to our 
results. Instead, we set up a test that checks whether the 
behavior we recover seems physically reasonable, and which 
allows us to look for any apparent numerical artifacts. Im¬ 
portantly, we want to check that the behavior that we see 
during the transition from an R-type to a D-type I-front is 
reasonable. Is the PDR in the proper position in relation 
to the ionization front, and does it evolves over time with 
the general expansion of the Hll region in the way that we 
would expect? While the ionization front remains R-type, 
the physical separation between the I-front and the edge of 
the PDR is very small and so we do not expect to resolve 
the PDR layer. However, once the I-front transitions to the 
D-type, we expect the PDR to broaden and so we should 
start to be able to resolve it. 

We use the parameters given in Table for the compo¬ 
sition of the medium and as input for the chemistry module, 
except for the initial gas and dust temperatures which we 
adjust to 10 K and 100 K, respectively, and fcuv, which is 
set from the properties of the radiation source. The simula¬ 
tion domain extends ±6 pc in each direction and we choose 
an initial number density of 1000 cm“®. Our setup corre¬ 
sponds roughly to the test case F2 in the previous section if 
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Figure 10. Profiles of the isothermal PDR tests F1-F4 after they converged in time. The light gray lines are the results from the eight 
PDR codes compared in |Rollig et al.||2007| l and the solid lines show the results obtained using Fervent. The surface densities are divided 
by the physical distance from the cloud edge to yield number densities. The density of the molecular cloud and the intensity of the 
incoming radiation field are annotated in the first row of the plot. The dashed lines show physical distances from the cloud edge. Our 
simplified chemistry is able to reproduce the results for the hydrogen species but only roughly approximates the carbon chemistry to an 
order of magnitude. 


we were to omit the radiation beyond the Lyman limit from 
the present test. 

We resolve the Stromgren sphere radius of 1 pc, which is 
only strictly defined for a purely atomic hydrogen medium, 
with a minimum cell size of 0.094 pc for the coarsest simula¬ 
tion with 128^ cells. This guarantees that the initial ionized 
region is captured so that the I-front expands properly. In 
addition, this high spatial resolution is needed to resolve the 
thin transition layer from molecular to atomic hydrogen. 

The radiation source is located at the center of the do¬ 
main and is fully characterized by its blackbody spectrum 
of Teft = 4 X 10"^ K and a stellar radius of 5 x 10^^ cm. 
The derived values for Afii, 2 , Nis.e etc. that are used in our 
radiative transfer scheme are shown in Table Finally, all 


of the physics described in Section is used in this test, 
including the effective cross-sections based on the effective 
blackbody temperature of the source, metal line cooling and 
the coupling of ionizing radiation to molecular hydrogen. 

In the remainder of this section, we look at the evolution 
of the combined ionization and photodissociation fronts at 
four times: during the R-type phase, at the transition time 
between R-type and D-type, during the early D-type phase, 
when the density contrast between the shell and the cavity is 
a factor of a few, and during the late D-type phase, when the 
density contrast between the shell and the cavity is large. 

Figurej^shows profiles of the temperature, density and 
pressure at these times. We find that the evolution of the 
combined fronts is independent of the spatial resolution. The 
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Figure 11. Profiles of the non-isothermal PDR tests V1-V4 after the gas has reached thermal and chemical equilibrium. The light gray 
lines are the results presented in |Rdllig et al.| l |2007| | and the solid lines show the results we obtain using Fervent. The surface densities 
are divided by the physical distance from the cloud edge to yield number densities. The density of the molecular cloud and the intensity 
of the incoming radiation field are annotated in the first row of the plot. The dashed lines show physical distances from the cloud edge. 
Our simplified chemistry is able to reproduce the results for the hydrogen species but only roughly approximates the carbon chemistry 
to within about an order of magnitude. The temperature we calculate lies within the spread of the other models. 
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Figure 13. Evolution of the thickness and chemical structure of the ionization and PDR fronts expanding in a uniform medium. We 
show the fronts at times of t = 1,100,500,1500 tree for spatial resolutions of 128^,256® and 512® cells from the top to bottom row. 
The black dots mark the position of the cell centers. The solid, dashed and dotted lines denote ionized, atomic and molecular hydrogen 
fractions of the gas. (Note that here we define xh 2 such that a value of one corresponds to fully molecular gas). 


Table 3. Radiation source parameters 


Symbol 

Derived value 

^11.2 

3.23 X 10"^® S-® 

iVl3.6 

1.61 X 10"^® s-i 

7715.2-1- 

4.70 X 10"^® s-i 

Ape 

1.51 X 10®® erg 

(Ell. 2 } 

12.33 eV 

ei3.6 

0.72 eV ® 

ei5.2+ 

2.94 eV 

eua 

4.44 eV 

0-13.6 

5.38 X 10"®® cm® 

015.2+ 

2.43 X 10“®® cm® 

OH 2 

6.01 X 10“®® cm® 


temperature of the Hll region is at around 7 — 8 x 10^ K, far 
below the value of 10 x 10^ K that is commonly assumed. At 
all times, the gas density contrast and thickness of the shell 
is limited by the spatial resolution. This leads to an ion¬ 
ized region that is slightly smaller for low resolution, as the 
shell always has a width of two to three cells, and the cells 
are larger in the low resolution simulations. In Figure |13[ 
we quantify the chemical evolution of the PDR and I-front 
structure by examining how the fractional abundances of 
the different chemical species change as we move across the 


structure. During the R-type phase, the ionized, atomic and 
molecular hydrogen layers are as thin as numerically possi¬ 
ble, i.e. only two to three cells wide. Gradually, the PDR 
spreads out which can be seen best in the thickness of the 
atomic hydrogen layer. At late times the dissociation region 
becomes several cells thick for all resolutions. However, due 
to the finite minimum cell sizes, they differ in total phys¬ 
ical size, i.e. the PDR structure is not yet converged. The 
transition of to CO follows closely the transition of H to 
H 2 , owing to our assumption that the CO photodissociation 
rate scales directly with the H 2 photodissociation rate. 


Finally, we also show slices through the simulation do¬ 
main in Figure for the highest resolution simulation with 
512® cells. We note that the gas in the PDR layer is heated 
to several hundred Kelvin. This leads to a reduced pressure 
contrast between the HII region and the ambient medium. 
The simulation was not set up to be in equilibrium, which 
results in a change in the temperature and pressure of the 
undisturbed gas over time. This gas cools from its initial 
temperature of 100 K to a temperature of a few tens of K 
at late times, close to its equilibrium value. 
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Figure 14. Slices of the combined PDR and ionization front expansion into a uniform medium. We show the density, pressure, temperature 
and atomic hydrogen fraction for four different output times. 


3.6 Photo-evaporation of a dense clump by two 
sources 

The environment in a typical star-forming region is not nni- 
form. It consists of large temperature and density contrasts 
with corresponding differences in the chemical structure. Af¬ 
ter a massive star forms in one of these regions inside a dense 
molecular core, it clears out its own vicinity of any remain¬ 
ing gas not used up in its creation. An initially isolated Hll 
region inside this core is formed which later breakes out and 
merges with others and affects the cloud as a whole. At this 
point, any remaining dense structures not host to the ini¬ 
tial star formation are impacted by the radiative feedback 
of several stars. 


We investigate the capability of our code to treat multi¬ 
ple sources illuminating a dense cloud in an idealized setup. 
The test consists of a uniform medium in a ±16 pc box with 
two identical sources at positions pi(x, y, z) = (—14, 0, 0) pc 
and P 2 {x, y, z) = (0, —14, 0) pc illuminating a homogeneous 
spherical over-density positioned at the origin with a radius 
of rdump = 4 pc. The parameters of the radiation sources 
are listed in Table We model two scenarios distinguished 
by the contrast between the clump density and the density 
of the ambient medium. In the first scenario, the hydro¬ 
gen nuclei number density in the ambient medium is set to 
n = lcm“^ and the temperature of the gas is set to 1000 K. 
The dense clump is initialized with a density a thousand 
times higher and a temperature of 10 K. 
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Figure 12. Expansion of a combined PDR and ionization front 
into a fully molecular uniform medium with a number density of 
1000 cm“®. The top panel shows the gas temperature and the 
small panel below it shows logarithmically scaled temperatures 
in the range of 10 to 1000 K to highlight the pre-heating due to 
the PDR. The third panel from the top shows density and the 
bottom pressure profiles from the source up to the position of the 
I-front for four different times, t = 1,100,500,1500 tree, which 
are shaded from light gray to black. The dotted, dashed and solid 
lines show results for spatial resolutions with 128^, 256^ and 512®, 
respectively. 


In these conditions, the initial recombination time in the 
ambient medium is ~ 0.12 Myr and each source creates an 
Hll region with a Stromgren radius of Rs = 58.3 pc which 
encompasses the whole simulation domain. The ionization 
fronts from both stars therefore reach the dense clump dur¬ 
ing their R-type phase, leaving no time for the ambient 
medium or the clump to react hydrodynamically. Accord¬ 
ingly, we only solve for the chemical and thermal structure 
within the first recombination time. This setup allows us to 
study whether there are any artifacts arising from our adap¬ 
tive ray-tracing scheme and rate calculation as well as the 
overlap and interaction of both ionization fronts. 

In the second scenario, we change the density contrast 
between the ambient medium and the clump to a factor of 
ten by increasing the density of the ambient medium density 


to n = 10 cm“® and decreasing the clump density to n = 
100 cm“®. This results in a recombination time of 0.012 Myr 
and the Stromgren radius decreases to Rs = 12.6 pc. These 
parameters allow us to study the hydrodynamical response 
of the ambient medium and the clump to heating by stellar 
radiation and the pressure from the surrounding ionized gas. 

Both setups have metal abundances and dust-to-gas ra¬ 
tios that are the same as in the PDR test (see Table 
and adiabatic indexes of 7 = 1.6667. We choose to initialize 
all gas to be fully molecular to highlight the effects of the 
PDR although realistically the ambient gas should be nearly 
atomic for both scenarios. 

Figure shows the evolution in the chemical compo¬ 
sition and the temperature of the gas for the first scenario 
for a series of different output times. The temperatures and 
chemical structure equilibrate after just 0.5 recombination 
times, with no expected change over the next few recom¬ 
bination times until the I-front reaches its D-type phase. 
During that phase, the temperature of the Hll region drops 
dramatically from around 28 x 10^ K in the ambient medium 
when the two ionization fronts meet to 7.5 x 10^ K at the 
quasi-steady state of the setup. 

The spike in the temperature at the midway point con¬ 
necting the two sources is the result of the added up heating 
and ionizing rates from both stars. Once both Hll regions 
fully merge and a large optically thin bubble is formed, the 
heating rate in this region becomes more similar to that 
near the sources and the temperature spike disappears. At 
all times the evolution of both I-fronts is symmetric, i.e. 
there is no biasing based on source position or the order in 
which the rays traverse the domain. 

At the surface of the dense clump, a PDR forms that is 
dominated by atomic hydrogen. Because of the high clump 
density, this PDR layer is narrow: the transition from fully 
ionized to fully molecular gas takes place in a distance of 
only around 1 pc. 

Very close to the radiation sources, the temperature 
drops by two thousand Kelvin at equilibrium compared to 
the temperature at intermediate distances. This artifact is a 
consequence of the fact that when solving the chemical rate 
equations, we require our ODE solver to produce accurate 
results only for those chemical species whose abundances 
exceed some specihed absolute tolerance. We do this on the 
grounds of computational efficiency - it makes little sense 
spending a large amount of computational time to accu¬ 
rately compute the Hi abundance if this is e.g. only 10“^® 
- but it means that when the ionization rate is very large, 
the atomic hydrogen abundance can become so small that 
it falls below this tolerance. If it does so, then the solver is 
at liberty to set it to zero. From the point of view of the 
chemical evolution, this introduces negligible error, as the 
gas is dominated by H'*'. However, it does affect the ther¬ 
mal evolution, since if xh is zero, there is no photoionization 
heating, whereas if xn is merely very small, the photoion¬ 
ization heating rate can be significant if the ionization rate 
is high. This artifact can be eliminated by reducing the ab¬ 
solute tolerance used in the solver, although this has the 
effect of increasing the computational cost of the entire sim¬ 
ulation. Alternatively, it can also be eliminated by putting 
a floor on the atomic hydrogen abundance, preventing it 
from ever dropping completely to zero, although again this 
has a measurable computational cost. However, in practice 
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Figure 15. High density contrast test of a dense clump irradiated by two sources. Shown is the evolution of the combined I- and PDR- 
front with no hydrodynamic evolution. The temperature profiles in the top row are taken along the dashed white line, the connecting 
line between both stellar sources. The profiles are plotted for resolutions of 128^ (dotted) 256^ (dashed) and 512^ (solid). The slices 
shown in the figure are taken from the highest resolution simulation. The profiles and slices are taken at three different times: when the 
I-fronts have just met, at an intermediate time when the ionization fronts have passed the clump, and at equilibrium. The spikes in the 
temperature profile at intermediate times are a transient feature introduced by the adaptively split rays, which introduce a small error 
in the sampling of the radiation field close to the source positions. 
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Figure 16. Errors and temperature structure after the high den¬ 
sity contrast setup reached equilibrium. The upper left plot quan¬ 
tifies the error from ray-splitting and the effect of spatial res¬ 
olution. The plotted temperature profiles are taken along the 
white dashed line for simulations with 128® (dashed green), 256® 
(dashed red) and 512® (dashed blue) cells. The relative error 
is computed from a control simulation with no adaptive ray¬ 
splitting at a resolution of 128®, denoted T^ef (solid black line). 
Iso-temperature contours of 7000, 7200 and 7400 K are over¬ 
plotted on the control run temperature slice in the top right panel. 
The bottom left panel shows the error incurred from adaptive ray¬ 
splitting in the temperature map and the panel on the bottom 
right shows a control run with no photo-electric heating. Note 
the different temperature scalings in the temperature slices. 


we only see this effect in extremely close proximity to very 
strong ionizing sources, and we do not find it to have any sig¬ 
nificant dynamical effects, suggesting that in more realistic 
simulations, this effect is likely to be harmless. 

Figure illustrates other numerical artifacts of our 
radiative transfer scheme. Some error is introduced by the 
adaptive splitting of rays that leads to a change in the sam¬ 
pling of the radiation field at the splitting radius. This re¬ 
sults in slight offsets in the accumulated gas columns and to¬ 
tal intersection length of all crossing rays which changes the 
calculated rates in comparison to an un-split case. We com¬ 
pare a control simulation with no adaptive splitting, where 
the initial number of rays is large enough to sample all cells 
sufficiently at all distances to one with ray-splitting. The er¬ 
rors shown in the upper and lower left panel in Figureare 
generally of the order of 1% at most. Another, transient fea¬ 
ture stemming from the adaptively split rays are the spikes 
observed in the temperature profiles at intermediate times 
in Figure |15| The amplitude of the spikes can be reduced 
by increasing the initial HEALPix level at additional compu¬ 
tational cost. For the default parameters the error around 
the source position is of 15 %. We tested the impact of the 
initial HEALPix level in control simulations utilizing the low 


density setup and found that it had no impact on the overall 
evolution. 


The fluctuations around the shadow cast by the dense 
clump and in the PDR stem from the varying number of rays 
that enter the cells. This depends on the current rotation 
of the HEALPix sphere, which leads to the calculated rates 
fluctuating slightly. As the edge of the shadow and PDR 
react very sensitively to changes in these rates, the error 
grows accordingly. 

Spatial resolution has an impact on the temperature at 
the source positions, since the smaller the volume of the cell 
with the source, the stronger the radiation field within the 
cell. Effectively, the grid is unable to resolve the point-like 
source which leads to a resolution dependence. 


The iso-temperature contours show that the heated re¬ 
gions close to the stars are not perfectly spherical. The origin 
of the lobe-like structures lies in the approximation used in 
modeling the photo-electric heating. It is a geometric arti¬ 
fact introduced in equation (52 I where the cell face area is 
used to recover a flux per area. To avoid it, one would have 
to calculate the area seen by the ray by using a projection 
onto the ray. We perform a control run with no photo-electric 
heating to make sure that there are no additional artifacts 
from other processes. As can be seen from the bottom right 
plot in Figure the resulting temperature field is nearly 
flat with no discernible lobes or other features. 


The evolution of the second scenario is shown in Fig¬ 
ure The ionization fronts meet after approximately ten 
recombination times, calculated for the ambient medium. At 
around the same time they transition to their D-type phase. 
Two shells around the Hll region can be made out in the 
density slice, one generated directly by the over-pressured 
ionized gas, the other from the heat deposited in the PDR, 
upstream of the I-front. The edge of the shell is unstable as 
can be seen along the cardinal directions where the outer, 
thinner shell breaks up. 

The shock front is reflected off the dense clump and 
after one hundred recombination times an interference pat¬ 
tern emerges from the incoming and the reflected shocks. 
The clump itself is symmetrically compressed from each side 
over the whole evolution, with a corresponding increase in 
temperature of a few thousand Kelvin on the surface, but 
no significant heating on its inside. 

In the low density ambient medium, a thick PDR forms 
which pre-heats the gas before the I-front reaches the same 
position. It can be easily made out in the atomic hydro¬ 
gen fraction slice shown in the bottom row of Figure |17| 
It becomes thinner as it moves into the over-dense region. 
Similarly, at later times the PDR preceding the I-front be¬ 
comes narrower as the shell becomes denser and the visual 
extinction increases over a shorter distance. 

Finally, Figure shows a spatial resolution study of 
the second scenario. As described in the previous section, 
the thickness and density of the shell depends on the spatial 
resolution. Additionally, this test shows that the stability of 
the shell is influenced as well: the smaller the cell size, the 
larger the density fluctuations. For all simulation runs, we 
find convergence in the position of the shells and fronts as 
well as in the temperature structure and the geometry of 
the compressed clump. 
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Figure 17. The low density contrast test case of the expansion of two Hll regions impinging on a dense clump with hydrodynamic 
response. Slices in density, temperature and atomic hydrogen fraction of the highest resolution runs with 512^ cells are shown. The 
evolution is plotted at three snapshots in time which correspond to recombination times of 1, 10 and 100 in the ambient medium. 


4 SUMMARY 

In this paper, we present Fervent, a radiative transfer code 
module for the magnetohydrodynamical adaptive mesh re¬ 
finement code FLASH 4. Fervent is designed to model the 
effects of radiation from massive stars without assuming ei¬ 
ther thermal or chemical equilibrium, as is otherwise com¬ 
mon in simulations of stellar feedback. It allows us to self- 
consistently evolve the chemical, thermal and density struc¬ 
ture of the interstellar medium surrounding a massive star 
over time. 


We are able capture the combined effect of the I-front 
and PDR front expansion necessary to accurately simulate 
Hll regions. This is in contrast to most dedicated photon- 
dominated region (PDR) codes that only calculate the final 
equilibrated state with no explicit time evolution and in a 
static density field. 


Fervent is based upon the ray-tracing scheme outlined 


in Wise & Abel (20111, which likewise utilizes the HEALPix 
library ( |G6rski et al.||2005| ) to adaptively split rays to sam¬ 
ple the radiation field. This approach is well-suited for com- 
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Figure 18. Resolution study for the low density contrast case. The late stage evolution after one hundred recombination times is shown 
for simulations with resolutions of 128^, 256® and 512®. 


billing with the adaptive mesh structure in FLASH 4, as de¬ 
scribed in detail in the Appendix. Our chemical and ther¬ 
mal treatment is designed for modelling the impact of stellar 
feedback on the present day ISM. We account for all of the 
main radiative heating processes: photoionization (of H and 
H 2 ), H 2 photodissociation, the vibrational pumping of H 2 
by FUV photons, and photoelectric heating by dust grains. 
We show that to properly capture these effects requires us 
to split up the photons into at least four frequency bins. 

We describe how to calculate photoionization and pho¬ 
todissociation rates in a fashion that is independent of spa¬ 
tial resolution, and show how to couple them to a fast bare- 
bones chemical network used to model the hydrogen chem¬ 
istry in the ISM. We also show how to treat a single chemical 
reaction that overlaps several energy bins and photons in a 
single energy bin that couple to multiple species in a photon 
conservative way. 

We test the code extensively from simplified to full se¬ 
tups where all physics modeling is included. Fervent repro¬ 
duces the well-known analytic results of I-front expansions 
in an atomic medium. In addition, we find very good agree¬ 
ment with results obtained from dedicated PDR codes for 
the structure of the hydrogen dissociation front as well as the 
temperature structure of this region. Disagreements arise for 
CO, which at present we can only approximate to within 
an order of magnitude owing to the simplifications made in 
our chemical network in order to decrease its computational 
cost. 

We explore under what conditions the hydrodynami- 
cal response to the thermal feedback is captured. We show 
that we have to limit the change in the hydrogen fraction 
to 10% over a single simulation time step and resolve the 
initial Stromgren-sphere by at least one resolution element 
for a proper expansion of the I-front. The effect of changing 
the spatial resolution is also thoroughly examined, and we 
show that apart from a few well-known issues to do with nu¬ 
merical diffusivity (shell thickness and stability, deviations 
from perfect spherical symmetry), all of our results converge 
in simulations with both low and high resolutions. 

In its current state Fervent is an extremely capable ra¬ 
diative transfer module with many potential applications in 
the field of star formation and the dynamics of the interstel¬ 


lar medium, which we plan to improve further in the future 
(e.g. proper treatment of the carbon chemistry). 
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APPENDIX A: RAY-TRACING 
IMPLEMENTATION DETAILS 

In the following, we describe the general structure of the 
adaptively split ray-tracing implemented in FLASH 4. 

A1 Data structures 

We distinguish two major data structures, the adaptive 
mesh used to decompose the computational domain and 
a two-dimensional array exclusive to the ray-tracing. The 
mesh structure is implemented as a fully-threaded tree that 
holds all of the hydrodynamical quantities for each resolu¬ 
tion element on each refinement level in the simulation. The 
ray array saves all information used in box-ray intersections, 
radiation source properties and some of the data gathered 
during sampling of the radiation held. Table [AT] shows the 
quantities carried in the mesh data structure that are ac¬ 
cessed and changed. Other quantities not relevant to the 
routine, e.g. pressure etc., are omitted. The data saved in 
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Table Al. quantities accessed and saved in 

name 

mesh data structure 

symbol 

gas density ^ 

p 

H abundance 

Xu 

H2 abundance 


abundance 


cell size 

d 

dissociation rate ^ 

fcuv 

ionization rate 

h- 

'^lon 

direct dissociation rate 

^dis 

ionization heating rate 

Fion 

direct dissociation heating rate 

Fdis 

photo-electric energy flux 

Fpe 


Table A3. Ray-tracing parameters 


Parameter name 

Default value 

pixel to cell face ratio (fa) 

4.0 

initial HEALpix level 

4 

maximum number of rays per comp, domain 

30 X 10^ 

rotate rays 

True 

split rays 

True 

stopping Ay 

10.0 

bundle rays 

5 

comm, interval 

20 

periodic in x 

False 

periodic in y 

False 

periodic in z 

False 

periodic box length 

1.00 


Table A2. quantities saved in ray data structure 
name 

direction, x component 

direction, y comp. 

direction, z comp. 

source position, x comp. 

source position, y comp. 

source position, z comp. 

number of ionizing photons in 

number of ionizing photons in £’15.2+ 

H ionizing photon energy in Eiz.Q 
H ionizing photon energy in -E15.2+ 

H2 ionizing photon energy 

average UV photon energy 

number of UV photons at source 

5.6 to 11.2 eV energy flux at source 

H ionization cross section for photons in Eiz.q 

H ionization cross section for photons in -E15.2+ 

H2 ionization cross section 

total H column 

H2 column 

HEALpix level 

HEALpix number 

total distance from source 

block identifier 

source identifier 

domain identifier 


symbol 

^dir 

ydir 

•^dir 

3^pos 

Vpos 

^pos 

^^15.2+ 

(Eis.e) 

(-E15.2+) 

(Sdis) 

(-^11,2) 

A^uv 

Fpe 

^ 13.6 


<^H2 

Nh 

Nh2 

Ih 

Dh 

rs 

bid 

®id 

did 


the ray array is shown in Table Besides these two main 
structures, auxiliary arrays are used in the parallel commu¬ 
nications routines, as described in Section |A4[ During the 
ray-tracing step, heating and reaction rates are saved into 
the tree structure for later use in the chemistry module. 
These quantities are shown in bold in Table [XT] 

The ray array is allocated at the start of the simulation 
and persists throughout the run. This has the advantage 
that we do not have to allocate and deallocate each array 
In each ray-tracing step, but has the disadvantage that we 
have to specify a maximum number of rays. 


A2 Algorithmic overview 

The ray tracing module is divided into several routines mo¬ 
tivated by the organization of resolution elements (cells) in 
uniform mesh blocks with extents Aceiis x Aceiis x Aceiia in 
the PARAMESH (MacNeice et al. 20001 tree data structure, 
the default AMR module in FLASH 4. Blocks, numbered by 
a locally unique identifier (ID), are the nodes of the tree and 


each time a block is marked to be refined, 2’" child blocks 
are created on the next higher level of the tree, where m 
Is the spatial dimension of the simulation. For blocks with 
2x2x2 cells this tree reduces to the commonly-used oct- 
tree. 

The main routine of the ray-tracing module Is called from 
Inside the chemistry module before the chemical network is 
solved, to make sure that the hydrodynamical data used in 
the chemistry and ray-tracing modules are the same. In prin¬ 
ciple, it could be called from any point during one timestep 
in the simulation. The main routine consists of two steps 
(see Figure [AT] |: the generation of rays with properties taken 
from a global source list, and a ray-tracing step, where the 
rays are advanced through the computational domains. 

For our ray-tracing module we assume a global list of all 
source positions and properties is available on all computa¬ 
tional domains. The source data structure carries the effec¬ 
tive temperatures and luminosities of the stars, from which 
the properties that describe the emitted radiation are calcu¬ 
lated (see Table A2|. Before any ray-tracing is performed, 
initial rays distributed on a HEALPix sphere are generated 
around each source location. The spheres are rotated by 
a randomly chosen angle every timestep to wash out any 
alignment artifacts with the Cartesian mesh structure. All 
generated rays are appended in the ray array without dis¬ 
tinguishing their origin. This has the consequence that we 
are not able to explicitly track individual contributions to 
the chemical state of a given cell from a given source. 

All domains containing a source begin ray-tracing the 
hrst ray in the ray data structure, while all others wait to re¬ 
ceive rays from neighboring domains. A ray is followed along 
through the blocks of the local domain it intersects until one 
of the following criteria is fulfilled: the ray leaves the global 
computational domain, it is stopped because the visual ex¬ 
tinction Av becomes large enough that there is no longer a 
need to follow the ray, or it reaches a block that is part of 
a neighboring domain and is queued for transport. Rays are 
moved between blocks by a tree lookup, where the proper 
neighbor is chosen during the intersection calculation. If the 
found neighbor block ID is on a different computational do¬ 
main, the corresponding ray is sent via message passing in¬ 
terface (MPI) communication to its target. Afterwards, the 
advancement routine moves to the next available ray, or if 
none are left locally, it stays in communication until all do¬ 
mains are done. 

The actual ray-box intersection calculation is done on a 
block level, where we take advantage of the uniform mesh 
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generate rays 


advance rays 



Figure Al. Structure of the ray-tracing module. The highlighted gray area marks the code for the traversal of a single ray. The routine 
name in bold is the only point in the module where outside data and data structures are requested and used. 


structure to optimize the ray-cell intersection (see next sec¬ 
tion and Figure |A2| for details). During the ray traversal 
of a single block, each time a new cell is entered the ray 
splitting criterion fn is evaluated. If the ray has to be split, 
the current slot in the ray array is overwritten by one of 
the child rays, and the additional three child rays are ap¬ 
pended to the end of the ray data structure. The splitting 
routine also checks whether any of the child rays move off 
the global domain, which might be the case if rays are split 
in the vicinity of the simulation domain boundaries. After 
splitting, the current cell iteration inside the block is ter¬ 
minated and a new iteration is called recursively with the 
updated child ray data instead of iterating to the next ray 
in the ray array. In this way, we always try to follow one 
ray or one of its child rays directly to the edge of the local 
domain, where the ray is sent without waiting for all local 
rays to traverse the local domain. 

The ray splitting, rotation of the HEALPix sphere and 
parallelization approach combine to randomize the order in 
which rays intersect cells in each timestep. Which and how 


many sources illuminate a cell, as well as the total num¬ 
ber of rays that enter it, is therefore unpredictable during 
the ray-tracing step itself. For multiple sources, to assure a 
properly converged result, one would have to take into ac¬ 
count all contributions from all sources per cell simultane¬ 
ously, call the chemical update during the ray-tracing, and 
re-distribute the attenuated photon fluxes onto the inter¬ 
secting rays, which are then moved to the next cell. 

Instead, as rays are indistinguishable between sources, 
we trust that the employed timestep criterion holds for mul¬ 
tiple sources. This can be seen if a cell illuminated by two 
identical sources on each side is considered. The rays enter 
and see the same chemical composition and calculate the 
same rates, which are added up to twice the value. Thus, the 
composition changes twice as fast and the timestep should 
adjust accordingly. A caveat is that we employ a corrector 
timestep, i.e. the timestep is reduced after a fast change was 
detected, and the first time two or more radiation fronts in¬ 
teract (or a radiation source turns on), the rate of change 
of the chemical composition will be underestimated, and 
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we may therefore employ a timestep that it larger than we 
would desire. As the validation study described in the main 
body of the text demonstrates, the errors arising from this 
are usually very small. Our ray-tracing approach does not 
suffer from any biasing based on the order that radiation 
sources are treated, simply because we calculate and sum 
all rates from all sources in one step. 

The ray-tracing module interfaces with the mesh data 
during the iteration over all intersected cells in a block via 
one single routine. It concentrates the physics modeling in 
the form of rate calculations in one place, and deposits all of 
the information required for the subsequent chemistry step 
in the appropriate fields in the tree data structure. This 
makes it easy to include additional physics, such as coupling 
to a more extensive chemical network, or adding additional 
energy bins. 

Finally, we give an overview of all parameters and their 
default settings that determine the behavior of the ray¬ 
tracing module (see Table A3l. The parameters can be 
roughly divided into two sets, one that describes how the 
radiation field is sampled and the other describing the com¬ 
munication pattern. The two most important parameters are 
fn, whose value approximately corresponds to the number 
of rays per cell and secondly the boolean flag that indicates 
whether rays are to be split at all. In some cases, especially 
for known geometries, it might be computationally more ef- 
hcient to hx the number of rays by specifying the initial 
HEALPix level, and not splitting the rays. In that case, one 
has to make sure that all cells are intersected by at least one 
ray. 


The parameters that control the behavior of the parallel 
communication are described in Section lAll Each simulation 
boundary can be set to be periodic, in which case a ray is 
allowed to leave a border and enter the opposite one. The 
ray is terminated if it traverses a total length given by the 
periodic box length parameter. 


A3 Box-ray intersection 

There are two stages to the intersection calculation of an 
inhnite line (ray) with an axis-aligned box (cell) in our 
ray-tracing module. First, once the ray enters a new block, 
the point the ray exits is computed by checking each block 
boundary individually. 



picking the set of boundaries that point away from the 
source, tnext = max(ti, tu), and finally choosing the shortest 
distance to one of the boundaries in the set, t = min(tnext). 
Here, we have denoted the source position by s, the upper 
and lower block boundary coordinates by b“ and b* and the 
direction vector by n. The exit point x is obtained by adding 
the scaled unit vector to the source position, x = t • n -|- s. 
Figure |A^ illustrates the block-ray intersection. 

Initially all rays start at the source with t — 0. Each 
time the advancement routine iterates to a new ray, the cur¬ 
rent position p of the ray at the beginning of the intersection 


block 



Figure A2. An example of a ray-block intersection. In this case 
the entry index is (3, 0) and the exit index is (0, 3). The short- 
dashed arrows indicate the order in which the cells are traversed. 


calculation is computed from the radius and source position 
saved in the ray data structure. 

Instead of using the source position in the intersection 
calculation, one could use an intermediate point by adding, 
e.g. the last intersection point in the previous block. How¬ 
ever, due to limits in numerical precision the distance cal¬ 
culated from s, Ts = t • n, typically does not match the one 
from adding up all previous block intersection lengths dr^ 
to the current block k, 




(A2) 


The small numerical drift is due to the square-root opera¬ 
tion necessary to obtain dr and can lead to the ray-tracing 
routine failing if both approaches are mixed. For exam¬ 
ple, depending on how close a ray comes to a vertex or 
edge, calculation of the ray position from might intersect 
the block while calculation from the summed intersection 
lengths might not. 

The computation of the exit point also allows us to look 
up the direction of the next block the ray enters during its 
traversal. The position of the smallest element in the tu or 
ti set gives the direction in the positive or negative x-, y-, 
or 2 -direction, respectively, and can be directly used in the 
tree lookup. 

The operation with the highest computational expense 
is the intersection length calculation for each cell inside one 
block. We optimize the operation by interpolating x and p 
onto the uniform mesh inside the block with cell size A®. 
This gives integer entry and exit indices, and each time a cell 
is crossed the index is incremented or decremented accord¬ 
ing to the cell face the ray leaves. The intersection length 


itself is computed from equation (All, where b is replaced 


with the cell boundary position c. We only consider upper 
or lower cell boundaries in each intersection calculation, de¬ 
pending on the direction of the ray. This halves the number 
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of computations necessary per cell. Additionally, we update 
and recompute t only for the boundary updated in the di¬ 
rection the ray left the current cell. We find the smallest of 
the new set tnext = and from this the intersection 

length dr = |min(tnext) • n|. We call the rates calculation 
routine each time dr is computed for a given cell. Once the 
ray enters the last cell of the block it intersects, the entry 
and exit indices become equal, and we stop the cell iteration. 
The integer stepping is highlighted in gray in Figure [A^ 
Every time the ray’s total distance property r^ is up¬ 
dated, the splitting criterion is evaluated. If the ray is split, 
one of the child rays is selected, and a new exit point is cal¬ 
culated. The block ray-trace is then called recursively and 
the ray-tracing continues. 


A4 Parallelization 


Arguably, the most performance relevant part of the ray¬ 
tracing module is its parallelization for use on distributed 
memory machines. Such machines consist of nodes, each con¬ 
taining several computing units, interconnected by a high 
bandwidth network. Each node has a fixed amount of mem¬ 
ory that holds part of the global simulation data, i.e. the 
local domain. The simulation data is usually decomposed 
by space, where each node is assigned some part of the total 
volume. 

Radiation from point sources is a difficult problem to 
parallelize, owing to the fact that sources are not evenly dis¬ 
tributed in space. A single local domain could potentially 
end up with most sources. This unbalances the work load, 
and in the worst case leads to a large fraction of the comput¬ 
ing units idling while they wait for the ones holding sources 
to finish their ray-tracing step. 

For our ray-tracing module, two parallelization schemes 
are possible. The first option is a synchronous communica¬ 
tion pattern, where we wait until all rays are traversed lo¬ 
cally before doing any communication. Rays that leave the 
local domain are saved in a buffer, and once all computing 
units finish their ray-tracing, a global communication step 
is performed. The buffered rays are sent and then another 
ray-tracing step is performed. This process is repeated until 
no rays are left in the global domain. For a single source, 
this can be seen as a domain-by-domain approach. First, 
the domain holding the source ray-traces, while all others 
wait. In the second step, the neighboring domains ray-trace, 
while again all others, including the source domain, wait. 
With this communication pattern, the work load moves out 
in shells around the source domain. From an implementation 
perspective, the module would be cleanly split into a ray¬ 
tracing step and a communication step, where all network 
sends are posted and completed. 

The second parallelization option is an asynchronous 
communication pattern. In this pattern, rays are bundled 
and sent off in intervals, without postponing communica¬ 
tion until all rays are traversed locally. The network sends 
are completed when the target computing unit calls the com¬ 
munication routine, usually after it has traversed a number 
of local rays itself, or no local rays are left. This is set by 


the communication interval parameter (see Table A31. In 
this communication pattern, the work imbalance is reduced 
by involving neighboring domains earlier in the ray-tracing. 
This leads to an overall decrease in computational time, as 


time is not wasted idling while waiting for the synchronous 
communication step. 

One additional important aspect is the amount of data 
sent during MPI communication. This sets the size of the 
communication buffers and in turn determines the number of 
messages that can be sent at the same time. We try to opti¬ 
mize this by only communicating the bare minimum amount 
of data necessary. Specihcally, the ray properties that char¬ 
acterize the ray direction, HEALpix level and unique pixel 
identiher, from which the direction vector can be looked up, 
the attenuated ionizing photon flux, the total and molecular 
hydrogen columns, the block and source identifier and the 
traversed distance of the ray from the source. The commu¬ 
nicated data is shown in italics in Table [A^ The other ray 
properties are locally reconstructed from the global source 
list. In total, we reduce the amount of information sent per 
ray to a third of what it otherwise would be. 
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